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Introduction 

Because of the interesting science which can be performed 
using a satellite attached by a very long tether to a mother 
vehicle in orbit, such as the Space shuttle, NASA will deploy 
TSS-1 (Tethered Satellite System) in 1992. A very long tether 
(20 km in this case) has the possibility of undergoing 
oscillations of several different types, or modes, and higher 
harmonics of these modes. The purpose of this document is to 
describe a method for detecting the amplitude, frequency, and 
phase (and predicting future motion in the steady state) of these 
modes, in particular, the skiprope mode, using tethered satellite 
dynamics measurements. Specifically the rotation rate data about 
two orthogonal axes, calculated from output from satellite 
gyroscopes, are used. The data of interest are the satellite 

pitch and roll rate measurements. 

Aside from understanding tether dynamics, one reason it is 
important to diagnose and predict the skiprope motion of the 
tether is related to satellite retrieval. The retrieval 
mechanism has a limited acceptance angle and the tether skiprope 
motion can cause angular excursions of the satellite beyond this 
angle. Several methods are available to damp the skiprope, but 
for these to be successful, it is necessary that the skiprope 
amplitude, frequency, and phase be known accurately enough and 
that future values of the parameters be predicted for the time of 

application of the corrective procedure. 

NASA has determined to use two methods to diagnose skiprope 
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properties and predict future values. One of these, a Fourier 
transform domain approach, is the subject of this notebook. The 
other is a time-domain, state-space method being developed by 
Martin-marietta. It is described elsewhere. Much of the 
development and testing of the Fourier algorithm and code has 
been done by the authors at the University of New Orleans. It is 
very important, however, to give full credit to the other 
contributors to the development effort. First and foremost is 
Mr. Stan Carroll of NASA Marshall Space Flight Center, whose help 
with this project was essential. Mr. Don Tomlin, Mr. Keith 
Mowery, Mr. Zack Galaboff, and others from MSFC have also 
contributed. From the University of Southern Mississippi we have 
received help from Dr. Grayson Rayborn and Dr. Sam Howard. 

The method which is described in this notebook can be 
modified to diagnose tether skiprope motion when the tethered 
satellite is spinning. This modification has been accomplished, 
but it is not discussed in this document. An addendum will be 
issued to describe it. There is also a version of the current 
code which contains a subroutine which gives the time history of 
the skiprope in the data observation window. That modification 
is also not described in this notebook. 

Included in this document are a section on Theory and 
Concept Definitions, in which some of the background theory and 
definitions for the method are discussed with references given 
for further reading. An algorithm outline follows, listing all 
the general steps of the program and including an overview flow 
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chart. For the main program and every subroutine, there are 
definitions of variables/ flags, descriptions of the input and 
output, and discussions of the code. The principal subroutine is 
WORK, within which the Fourier transform and related calculations 
are performed. Other subroutines which are either called by WORK 
or the main program are HANN, which applies the HANN window; 

MEAN, for mean removal; F0UR1, for the calculation of the fast 
Fourier transform; LSCF, which performs a least squares curve 
fit; FBAND, which defines the search band in the transform domain 
to find the maximum of the transform magnitude; and READINDATA 
and ODTF for input/output . Also included in the section on 
program description is a quick reference which gives required 
external files, operator inputs, and a listing of warnings/abort 
situations. Following the program description is a section which 
gives program development history. This section is written both 
to show the trials which led to the selected algorithms and to 
inform the reader about other methods which did not perform as 
well as those adopted. Finally, in the main body of this 
document is a complete description of the test plan which was 
executed at the University of New Orleans. 

There are two Appendices. The first lists the program code 
which accomplishes skiprope observation for slowly varying 
skiprope parameters and prediction of future parameter values for 
steady state tether motion. Appendix 2, which contains the test 
code and results, forms the bulk of this document. It is so long 
that it has its own table of contents. It is broken up into six 
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parts, Appendices 2. A through 2.F, Appendix 2. A contains the 
programs for generating model signals. In Appendix 2.B are the 
results of the ECR verification table. Appendix 2.C lists the 
programs for ECR testing, while Appendix 2.D lists the programs 
for systematic testing. Appendix 2.E gives the results of 
simulation tests, and Appendix 2.F gives the systematic test 
results at Station 2 and Station 1. 
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Frequency Domain Skiprope Observer 
Algorithm Outline 

1) Access the data buffers produced by the preprocessor 
containing the pitch and roll gyro data. 

2) Select a time window of data to be analyzed. 

For each of the axes perform the following: 

3) Apply Hann window to the data. 

4) Calculate the mean. 

5) Remove the mean from the original data. 

6) Apply Hann window to the mean-removed data to reduce 
sidelobes in the Fourier transform. 

7) Use the FFT to calculate the real and imaginary parts of the 
Fourier transform of the Hann-windowed data. 

8) Calculate the amplitude spectrum of the transform. 
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Frequency Domain Skiprope Observer 


Algorithm Outline - continued 

9) Search the predefined frequency band to find the maximum in 
the amplitude spectrum. 

10) Fit a quadratic to the seven points centered around and 
including this maximum. 

11) Use this quadratic to define more accurately the largest 
value in the amplitude spectrum and its corresponding frequency. 

12) Calculate the satellite motion amplitude from this maximum 
amplitude value by a simple conversion. The frequency is the 
frequency of that motion. 

13) Use quadratic interpolation to find the real and imaginary 
values of the Fourier transform at the above calculated 
frequency. 

14) Use real and imaginary parts to define the phase of the 
satellite motion based on the model cos (2 n f t + <p) . 

15) Calculate the time index for the first maximum of the gyro 
signal w within the data window. 
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Algorithm Outline - continued 
CALCULATIONS UPON RETURN FROM WORK SUBROUTINE 

16) Average the x and y frequencies - call FAVG. 

17) Find the reciprocal of FAVG (average period) , call TAVG. 

18) Calculate constant WK = TETHER LENGTH * TAVG / (360 * PI) 

19) Calculate polarity: 

a) Compare time indices of maximum. 

b) For x time greater than y time. If time difference is 
less than 1/2 period, polarity is positive, else polarity is 
negative . 

c) For y time greater than x time. If time difference is 
less than 1/2 period, polarity is negative, else polarity is 
positive . 

20) For both x and y axes, calculate omega values: 

W (X or Y) = AMP (Y or X) * COS (2.0 * PI * FAVG * DT * (I - 1) 

+ PHASE (X or Y) ) 

21) For both x and y axes, calculate motion amplitude values 
U (or V) = - POLARITY * WK * W (X or Y) 

22) Calculate time values TIME = TO + (START INDEX +1-2) * DT 
YAW MANEUVER CALCULATIONS 

24) Time of maximum x value: TXMAX = DT * (INDEX OF MAX - 1) 

25) Time of maximum y value: TYMAX = TXMAX + 0.25 * TAVG * 

TSHIFT + TO 

(Note: TXMAX is a relative time, TYMAX is a mission elapsed 
time . ) 

26) Time for orbiter maneuver (burn time) : TT = TYMAX + TAVG * 
(INTEGER) 

(Burn times are integral numbers of periods from TYMAX.) 
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PROGRAM DESCRIPTION 


Main Program 


Definition of Variables/Flags 


Parameters 

NDIM - dimension of arrays 
NFT ~ number of points in 


X, Y, and TLG (set to 3000) 
the Fourier transform (set to 


HPNT -dumber of points used in the least squares curve 
fitting polynomial (set to 7) 


Constants 
PI - value of pi 
DFR - conversion 
57.2957795) 
RFD - conversion 


(set to 3.1415926) 

factor from radians to degrees 

factor from degrees to radians 


0.017453292) 


(set to 
(set to 


REPL? C - character variable of length 1 (replies by the 
onerator to prompts, either y, Y, n, or N) 

POLARITY ?cha?acter P variable of length 8 (either pos- 

itive or necrative) 


Variables read in from external file PARAM.DAT 
PF SRCH BAND - percentage used in calculating the freq- 
“ uency search band) # 

R ARM “ distance from orbiter C*M. to the center me 
of the deployer boom . . 

ODF TIME - logical flag to control printing of the time 

- data to output file F_TXYUV.DAT (initially 
set to .FALSE., i.e., don't print) 

ODF FFT - logical flag to control printing of frequency 

- domain data (FFT's) to the output files 
F_FFTX.DAT and F_FFTY . DAT (initially set to 

DENSITY - tether ' density^in^kg/km (set to 8 35 kg/km, 

--SSitTSS W &? tTsSiS r eters) 

MORB - orbiter mass in kg (set to 100, 00i o°km? 

ALTKM - orbit altitude in km (set to 325.0 km) 


Variables read from telemetry preprocessor (or external 

file IDFTXY.DAT) . 

TO - time tag for first point in buffer 
TF - time tag for last point in buffer 
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X(NDIM) - X axis gyro data array (deg/sec) 

Y(NDIM) - y axis gyro data array (deg/sec) 

TLG(NDIM) - tether length array (kilometers) 
mode - integer variable, signifying the ancsmode for 
the last time point 'TF' 
amcsmode = 0 indicates no valid data 
amcsmode = 1 indicates passive case 
amcsmode = 2 indicates yaw hold 
amcsmode = 3 indicates spin case 
LF - number of time points stored in each array 
M FLAG - logical flag to denote if the amcsmode 
" chlnged from 1 or 2 to a 0 or 3 between the 

times of TO and TF 

AMCSMOD^- 3 integer variable set to the value of MODE 
- starting time index (requested of operator) 

LE - last time index (requested of operator) _ 

LEB - total number of data points processed (LEB LE 

LB +1) 

DT - average sample time = (TF 
DF - frequency sampling — 1.0/(NFT uij 
FLOW - low frequency of the frequency search band 
FHIGH - high frequency of the frequency search ban 
S? - estimateof the number of data points to use 
for a minimum of 3 cycles of skiprope 
LEAST = INT(6.0/ (FLOW + FHIGH)) 
ttmcth — tether length used in calculations, firs 
TLHGTH estimated Iron 0.5* (TLG ( 1) t TLG(LF) ) , and 
finalized as 0 . 5 * (TLG (LB) + TLG(LE)) 
tshift - time shift from first point in buffer, i.e., 

TO, and the first time point used in the run 

TSHIFT = DT*(LB - 1) . , 

TMIDPT - time point of the middle of the data win ow 
TMIDPT = TO + DT* ( (LB + LE -2)/2) 
e FLAG - logical flag to control yaw maneuver calcula 
S_FLAG logi^ ( o n ly*set to .TRUE, by direct operator 

reply of 'y' or 'Y* after prompt) 

PSIGN - numerical sign of the polarity (calculated) 

PSIGN = 1.0 indicates positive polarity 
PSIGN = -1.0 indicates negative polarity 
PSIGN = 0.0 indicates inability to predict the 
skiprope frequency - 

FREQX - calculated value of the skiprope frequency from 
the x axis gyro data „ opnY 

AMPX - calculated value of the x gyro rate at FREQX 
PHASEX - phase of FREQX relative to LB 
TWXMAX - time index of the maximum x gyro rate value 
52“*? me an™ value of Hann-windowed x axis gyro rate data 
g ?Lgx ! logical flag indicating whether the x axis 
~ values are good, i.e., G_FLAGX is see ro 

.TRUE, if FREQX is within the specified fre- 
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quency search band 

FREQY — calculated value of the skiprope frequency from 
the y axis gyro data 

AMPY - calculated value of the y gyro rate at FREQY 
PHASEY - phase of FREQY relative to LB 
IWYMAX - time index of the maximum y gyro rate value 
AVGY - mean value of Hann-windowed y axis gyro rate data 
G FLAGY - logical flag indicating whether the y axis 
— values are good, i.e., G_FLAGY is set to 

.TRUE, if FREQY is within the specified fre- 
quency search band 

FAVG - average skiprope frequency = 0.5*(FREQX + FREQY) 
TAVG - period of the average frequency = 1.0/FAVG 
WK - conversion factor from gyro rate values to ampli- 
tudes WK = 1000 . 0 *TLNGTH*TAVG/ (360*PI) 

UMAX - maximum in plane skiprope amplitude 
UMAX = WK*AMPY 

VMAX - maximum out of plane skiprope amplitude 
VMAX = WK*AMPX . 

TEST — time required to move from the x axis to the y 
3 l x is ( i n s 6 c J 

If IWYMAX .GT. IWXMAX, TEST = DT* (IWYMAX-IWXMAX) 
If IWXMAX .GT. IWYMAX, TEST = DT* ( IWXMAX- IWYMAX) 
TWX - time domain skiprope signal for the x axis (with 
out proper amplitude scaling) 

TWX = COS (2.0 * PI * FAVG * DT * (1-1) + PHASEX) 
TWY - time domain skiprope signal for the y axis (with- 
out proper amplitude scaling) 

rpyiY — COS (2.0 * PI * FAVG * DT * (I“l) + PHASEY) 
wx - time domain skiprope signal for the x axis (with 
proper amplitude scaling) WX = TWX * AMPX 
WY - time domain skiprope signal for the y axis (with 
proper amplitude scaling) WY = TWY * AMPY 
U — in plane skiprope amplitude (in meters) 

U = -PSIGN * WK * AMPY * TWX 
V - out of plane skiprope amplitude (in meters) 

V = -PSIGN * WK * AMPX * TWY 
time tag = TO + (LB + I - 2) * DT 


T - 
RNROT - 


TXMAX - 


TYMAX - 


number of rotations the orbiter should execute 
RNROT = ( AMIN1 (UMAX, VMAX) ) / (2.0 * R_ARM) 
time of maximum x axis gyro rate (should corre 
spond to when the tether is over the orbiter 
nose) = DT * (IWXMAX - 1) 

time of maximum y axis gyro rate (should corre 
spond to when the tether is over an orbiter 
wing) = TXMAX + 0.25 * TAVG + TSHIFT + TO 
TT - predicted times to execute the yaw maneuver 

TT = TYMAX + (K-l) * TAVG (K is an integer that 
runs from 1 to 5) 

OYAWANG — orbiter yaw axis angle in degrees 

BTDEL - burn time delay = PSIGN * OYAWANG/ (360 . 0*FAVG) 
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Common Variables 

LB, LE, DT, DF, FLOW, FHIGH, 

MORB, ALTKM 


DENSITY, TOTALL, MSAT , 


Input /Output Files 


All input and output files are opened and closed in the 
main program. 

Input Files (External) 

PARAM.DAT (read in the main program) „ 

The external input file PARAM.DAT consists of two lines 

having the following format: 

PF SRCH BAND , R_ARM, ODF_TIME, ODF_FFT 
DENSITY, TOTALL, MSAT, MORB, ALTKM 

IDFTXY.DAT (read in the subroutine READINDATA) 

The external input file IDFTXY.DAT simulates the pre- 
processed telemetry data stream. It consists of a 
maximum of 3000 lines with the following format: 

TIME, X(I), Y (I) , TLG(I) , MODE 

Output Files , 

F FFTX.DAT (written in the subroutine WORK) 

The FFT of the x axis gyro rate data is written to file 

F FFTX • DAT * 

F~FFTY.DAT (written in the subroutine WORK) . 

The FFT of the y axis gyro rate data is written to file 

F FFTY DAT « 

Both F^FFTX.DAT and F_FFTY.DAT have the following 

FREQUENCY, MODULUS, REAL PART, IMAGINARY PART 

Writing to both F_FFTX.DAT and F_FFTY.DAT is controlled 

by the logical flag ODF_FFT. 

F TXYUV.DAT (written in the main program) 

Time domain data is written to file F_TXVUV.DAT. The 
logical flag ODF_TIME controls writing to F_TXYUV.DAT. 
The format for F_TXYUV.DAT is: 

T, WX, WY, U, V (defined above in variable list) . 

F YAWMAN.DAT (written in the main program) 

If yaw maneuver calculations are done, the resuits are 
written to the file F_YAWMAN.DAT, which has the format. 
TMIDPT TF 

(K-l) /POLARITY, 360.0 * FAVG, RNROT, TT 

where K runs from 1 to 5 (other variables defined as 

above in the variable list) . 
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F_RECORD.DAT (written in subroutine ODTF) 

See the subroutine ODTF for a description of this file. 


Subroutine Calls (in order of calling) 


READINDATA ( TO , TF , X , Y , TLG , JMODE , LF , M_FLAG) 

FBAND ( FLOW , FHIGH , TLNGTH , PF_SRCH_BAND) 

Subroutine FBAND is called twice, the first time to 
return values of FLOW and FHIGH to use in estimating 
the number of data points necessary for 3 cycles of the 
skiprope, and the second to actually calculate FLOW and 
FHIGH for the frequency search band. 

WORK (11, X, AMPX, PHASEX, FREQX , IWXMAX, G_FLAGX, AVGX, 
ODF_FFT) 

WORK (12, Y, AMPY, PHASEY, FREQY, IWYMAX, G_FLAGY, AVGY , 
ODF_FFT) 

ODTF (TO, TMIDPT, TF, DT, LE, LB, LEB, TLNGTH, AMPX, 
FREQX, PHASEX, AMPY, FREQY, PHASEY, FAVG, TAVG, 

WK, PSIGN, UMAX, VMAX, FLOW, FHIGH, AVGX, AVGY) 

Code Discussion 


The parameters NDIM, NFT, and NPNT are set to the 
values 3000, 8192, and 7, respectively, the constants 
PI, DFR, and RFD are set to 3.1415926, 57.2957795, and 
0.017453292, respectively, and the arrays X, Y, and TLG 
are dimensioned to NDIM. AMCSMODE is declared as an 
integer, REPLY and POLARITY as character variables, 

G FLAGX, G_FLAGY , S_FLAG , M_FLAG, ODF_TIME, and ODF_FFT 
as logical variables, and PF_SRCH_BAND, R_ARM, DENSITY, 
TOTALL, MSAR, MORB, and ALTKM as reals. The variables 
LB, LE, DT, DF, FLOW, FHIGH, DENSITY, TOTALL, MSAT, 
MORB, and ALTKM are declared common. 

The external file PARAM.DAT is opened as logical 
unit 10. The values of PF_SRCH_BAND, R_ARM, 0DF_TIME, 
ODF_FFT , DENSITY, TOTALL, MSAT, MORB, and ALTKM are 
read and PARAM.DAT is closed. All output files are 
opened, with F_FFTX.DAT as unit 11, F_FFTY.DAT as unit 
12, F_TXYUV.DAT as unit 13, F_YAWMAN.DAT as unit 17, 
and F_RECORD.DAT as unit 18. 

The subroutine READINDATA is called to read in the 
preprocessed telemetry data. (At present, this data is 
simulated in the file IDFTXY.DAT.) READINDATA returns 
the values of TO, TF, JMODE, LF, M_FLAG, and the arrays 
X, Y, and TLG. The variable AMCSMODE is set to the 
value of JMODE. If M_FLAG is .FALSE., the operator is 
warned that the AMCSMODE changed during the data stream 
and the number of data points may be reduced. If the 
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value of AMCSMODE is either a 0 (indicating an invalid 
data set) or a 3 (spin case), the operator is alerted 
and the program aborts. 

The yaw maneuver flag S_FLAG is set to .FALSE. The 
operator is asked whether yaw maneuver calculations are 
required or not, and is advised that the calculations 
w ill not be performed unless requested. Only if the 
operator replies with ' y ' or 'Y' will S_FLAG be set to 
.TRUE, and yaw maneuver calculations executed. 

The average sample time DT = (TF - TO)/(LF - 1) 
and a preliminary tether length TLNGTH = 0.5 * (TLG(1)+ 
TLG(LF)) are calculated. The subroutine FBAND is 
called (passing this TLNGTH and P F_SRCH_B AND ) to return 
values of FLOW and FHIGH used in estimating the number 
of data points necessary to comprise 3 skiprope cycles, 
LEAST = INT(6.0/ (FLOW + FHIGH)). LEAST is printed to 
the screen, and the operator is prompted to enter LB, 
the starting time index, and LE, the last time index. 

If LE - LB + 1 is an even number, set LE = LE - 1 so 
that the total number of time points LEB = LE - LB + 1 
is odd. (An odd number is necessary for proper use of 
the HANN window subroutine called by the WORK subrou- 
tine ) The tether length TLNGTH is recalculated as 
TLNGTH = 0.5 * (TLG(LB) + TLG(LE)). This value is 
printed to the screen for the operator's approval. 
Subroutine FBAND is called again with the new value of 
TLNGTH and returns the values of FLOW and FHIGH used as 
the end points of the frequency search band. These 
values of FLOW and FHIGH are printed to the screen for 
the operator's approval. (The operator may change the 
values of TLNGTH, FLOW, and FHIGH if disapproved.) The 
values of LB, LE, LEB, and TLNGTH are printed to the 
screen. DF = 1.0/(NFT*DT) (the frequency sample rate), 
TSHIFT = DT * (LB - 1) (the time shift from TO, the 
first point in the data buffer, to the time of LB, the 
start index), and TMIDPT = TO + DT * ( (LB + LE - 2)/2) 
(the time of the midpoint of the data window) are now 
calculated. 

The WORK subroutine is called twice, once passing 
the array X and the flag 0DF_FFT , and once passing the 
array Y and the flag ODF_FFT. WORK returns the values 
of AMPX, PHASEX, FREQX, IWXMAX, G_FLAGX, and AVGX after 
the first call, and AMPY, PHASEY, FREQY , IWYMAX, 

G FLAGY, and AVGY after the second call. PSIGN is set 
to the default value of 0.0 (if PSIGN remains as 0.0 
then this indicates failure to predict the skiprope 
frequency) . The flags G_FLAGX and G_FLAGY are checked, 
with 4 resulting cases: 

1) If both G FLAGX and G_FLAGY are true, calculate FAVG 
as the avirage of FREQX and FREQY, the period TAVG 
as the reciprocal of FAVG, the constant WK = 1000.0* 
TLNGTH*TAVG/ (3 60*PI) , UMAX = WK*AMPY, and VMAX — WK* 
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AMPX. Print the values of AMPX, PHASEX, and FREQX, 

AMPY, PHASEY, and FREQY, 57.3 * (PHASEX - PHASEY) 

(the phase difference between x and y) , UMAX and 
VMAX , and TAVG to the screen. . = 

if G FLAGX is true and G FLAGY is false, set FAVG 

’ FREQx! TAVG = 1.0/FAVG, WK = 1000.0 * TLNGTH * TAVG / 
(360*PI) VMAX = WK*AMPX, and UMAX = 7777.0. Print 
to the screen the warning that the y axis data is 
suspect, with the calculated frequency outside the 
search band. (The frequency returned for y is the 
predicted midpoint of the search band.) This data 
should not be used without caution. Neither the 
polarity nor the yaw maneuver calculations are per- 
formed. The value of VMAX is printed to the screen. 

Tf G FLAGY is true and G FLAGX is false, set FAVG 
FREOY TAVG = 1.0/FAVG, WK = 1000.0 * TLNGTH * TAVG/ 
(360*PI) , UMAX = WK*AMPY , and VMAX = 7777.0. Print 
to the screen the warning that the x axis data is 
suspect, with the calculated frequency outside the 
search band. (The frequency returned for x is the 
predicted midpoint of the search band.) This data 
should not be used without caution. Neither the 
polarity nor the yaw maneuver calculations are per- 
formed. The value of UMAX is printed to the screen. 

4 \ if both G FLAGX and G_FLAGY are false, set UMAX, 

VMAX, FAVG, TAVG, WK, and PSIGN to 7777.0 and print 
to the screen that both axes are bad and offer 
suggestions for action: 1) look at the time plots o 
the gyro signals; 2) look at the FFT plots; an ) 

widen the search band. ^ . . .. 

The polarity PSIGN is calculated using the dif 
ference between IWXMAX and IWYMAX and competing to 
0 5* TAVG. Two cases are checked: 1) IWYMAX greater than 
IWXMAX, and 2) IWXMAX greater than IWYMAX. For case 1) 
TEST is set to DT * (IWYMAX - IWXMAX). If TEST is 
areater than 0.5*TAVG, then PSIGN = -1.0, else PSIGN - 
1.0. For case 2) TEST is set to DT * (IWXMAX - IWYMAX). 

If TEST is greater than 0.5*TAVG, then PSIGN — 1.0, else 

PSIGN = -1.0. Print PSIGN to the screen. 

If the logical flag ODFJTIME is set to .TRUE., then 
calculate TWX, TWY, WX, WY, U, V, and T, and print T, 

WX WY, U, and V to the file F_TXYUV.DAT. Note that all 
S' theie variables are only calculated if ODFJTIME. is 
true. (As stated in the variables definition section, 
TWX = COS (2.0 * PI * FAVG * DT * (I - 1) + PHASEX), 

TWY = COS (2.0 * PI * FAVG * DT * (I - 1) + PHASEY), 

WX = TWX* AMPX, WY = TWY* AMPY, U = -PSIGN*WK*AMPY*TWX, ^ 

V = -PSIGN*WK*AMPX*TWY , and T = TO + (LB + I - 2) DT.) 

If the yaw maneuver logical flag S_FLAG is .TRUE., 
then calculate the number of rotations the orbiter 
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should execute, RNROT = (the minimum of (UMAX, VMAX) ) / 
n * R arm). If PSIGN = 1.0, set the character 
variable - POLARITY to 'POSITIVE', and if PSIGN = -1.0, 
set POLARITY to 'NEGATIVE'. Calculate the time when 
the tether is over the orbiter nose, TXMAX D 
( IWXMAX - 1) (note that this time is a relative time 
£*? beginning of the data window only!) . 9*^^ 

the time when the tether is over an ° rb i be ? 

= TXMAX + 0.25 * TAVG + TSHIFT + TO (note that this is 

an absolute or mission time quantity) . ^«t th 

values of the data window midpoint time TMIDPT and tne 
last time in the data buffer TF to both the screen and 
the file F YAWMAN.DAT. The time TYMAX must be adjusted 
to account-for the orbiter orientation with respect to 

verify'^the'orblter £?££ an g ?e°S?AW»G° Osculate the 
burHime delay btdL - PSIGN * OOHC/ (360.o*favg) and 
tvmax Compare the time TYMAX to Tr • 

M TYMA^is* less^than or equal to TF. add multiples of 

^r^? d eer?F?° ITi S. T c“cutate V ?L U ?iL 

revolution label (K - 1) , POLARITY, 360.0 * FAVG, 

RNR ° T Regardless of whether the file F TXYUV.OAT has been 
printed or not, irrespective of the value of S FLAG, 
for all 4 cases of G FLAGX and G FLAGY combinations, 
call sibroS?!ne 0 DTF"and write the fiLe F RECORD.DAT. 

TLNS?H th LlPX 1U FREQX; PHASEX^^Y, FREQY , PHASEY, FAVG, 
TAVG"psiG™ E Skx, VMAX , FLOW , FHIGH AVGX AVGY. 
Close all files (logical units 11, 12, 12 , 1 , 


Subroutine WORK(IOA, ANG, AMP, PHASE, FREQ, ITMAX, 

G_FLAG, BIAS, FFT_FLAG) 


Definition of Variables/Flags 


Variables passed as arguments . ffl 

IOA - logical unit number for writing output 

unitll is for file F_FFTX.DAT, 12 for F_FFTY.DAT 
ANG — ayro rate data (either x or y axis) , , 

AMP - amplitude of the gyro rate data at the calculate 

skiprope frequency FREQ 

PHASE - phase of the skiprope frequency relative to the 
beginning of the data window at index LB 
FREQ - calculated skiprope frequency 
ITMAX - time index of the maximum gyro rate 
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G FLAG - logical flag indicating . wheth ®^^® r ®REQ n is 
value of FREQ is good, i.e., whether FREQ l 

found within the bounds of the frequency 

BIAS - mea^valuiTof the gyro rate data array ANG after 

EFT fla^- 1 logical flag controlling the writing of the 
— files F FFTX.DAT and F_FFTY.DAT; initially 

set to • TRUE • 

NDIM^dimension of the arrays ang and aux (set = 3000) 
UrnTM - dimension of the complex array awo (set - 8200) 
S““ numbe?of points in the Fourier Transform (set = 

NPNT - 8 number of points used in the least squares curve 
fitting polynomial (set — 7) 

Constants 

DFR"- S Sonversion 1 factor from radians to degrees (set to 

RFD - conversion 5 factor from degrees to radians (set to 
0.017453292) 

Common iVariables index (of gyro rate data array ang) 

LE - last time index (of gyro rate data array ang). 

DT - average sample time 

fLw £ - r !r?req;Jn?y e bo2n? of the frequency search band 
FHIGH - high frequency bound of the frequency searc 
band 

Other Variables arrav shifted so first index is 1 

ACX - » r ° rfe data array^ Shit Fourier Transform 

AWO - complex data array, useu LUJ - _ _ 

NTB1 - number of data points in array ( . 

rii - index shift used in creating array AUX (1 LB) 

- maximum value of gyro rate data found in array 

IFRST - A i£dex of the transformed array AWO corresponding 
to FLOW (IFRST = 1 + INT (FLOW/DF) ) 

I LAST - index of the transformed array AWO corresponding 
to FHIGH (ILAST = 1 + INT (FHIGH /DF) ) 
fr - frequency corresponding to index I in transformed 
array AWO (FR — (I 1) *DF) 

Vi? — index of maximum modulus of array AWO , _ 

XFREQ 1 - array Sf 7 points, consisting of the moduli of 
XFRE0 array of the array Aw0 with indices^ 

centered about KF, used in the least squares 
curve fitting of FREQ 
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PHIMAG 


PHREAL 


FQ_PO - 

PHASEI 
PHASER 
SCALE • 


_ array of 7 points, consisting of the imaginary 
parts of the 7 entries, of the array AWO with 
indices centered about KF, used in the least 
squares fitting of PHASEI 

- array of 7 points, consisting of the real 
parts of the 7 entries of the array AWO with 
indices centered about KF, used in the least 
squares fitting of PHASER 

interpolated index value returned by the curve 
fitting subroutine LSCF 

- imaginary part of PHASE, interpolated by LSCF 

- real part of PHASE, interpolated by LSCF 
scaling factor to give transformed data m units 
of deg/ sec and represent actual rate data 


Output Files 


F FFTX.DAT, F_FFTY.DAT _ . . „ , 

— File F FFTX.DAT has logical unit number 11 (storea 
in IOA) and F FFTY.DAT has logical unit number 12. For 
long tether leHgths (small skiprope frequencies) , 1006 

lines are printed; for short tether lengths (larger 
skiprope frequencies) , 336 lines are printed. Each line 


has the format: 

FR, XMAX, REAL ( AWO ( I ) ) , AIMAG (AWO (I) ) 
(Here XMAX = SCALE * modulus of AWO(I)) 


Subroutines Calls (in order of calling) 


HANN (NTB1, AUX, BIAS) 

F0UR1 (AWO, NFT , 1) 

LSCF (FQ P0, XMAX, XFREQ , 1, G_FLAG) . 

G FLAG il s4t after this first call to LSCF (with the 1 
in the 4th argument) . If G_FLAG is .TRUE. , then LSCF is 
called twice more to interpolate the imaginary and real 


parts of PHASE: 

LSCF (FQ P0, PHASEI, PHIMAG, 2, G_FLAG) 
LSCF (FQ_P0, PHASER, PHREAL, 2, G_FLAG) 


Code Discussion 


The values of the gyro rate data array ANG, logical 
unit indicator IOA, and file print flag FFT_FLAG are 
passed in the calling statement, as are the common van 
ables LB, LE, DT, DF, FLOW, and FHIGH. The parameters 
NDIM, NCDIM, NFT, and NPNT are set to the vaiues 3000, 
8200, 8192, and 7, respectively, and the constants PI, 
DFR, and RFD to 3.1415926, 57.2957795, and 0.017453292, 
respectively. The real array AUX is dimensioned to 
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are declared logical variables. 

The number of data points NTB1 in the array ANG 
fNTBl - LE - LB + 1) and the first index shift LB1 (LB1 
i i - LB) are calculated. Letting the index I range 

LB Senate th4 mean of the array AUX (and thus of 
ANG) and apply the HANN window to AUX by calling the 

SUbr °Make e tht N complex array AWO by setting the real 
marts of the first NTB1 entries of AWO equal to the 

corresponding entry in AUX, with the imaginaryjarts 

S pfnd t ?he-°Fiu^er P ? d ran2?orm e of e ^0°using the FFT routine 

F ° UR1 Search i f or S the^maximum ,, modulus a of R the P transformed 

set XMAX = 0.0, and calculate the 
array AWO. F ^t, set; - 1 + INT (FLOW/DF) and I LAST = 
frequency indices IFRST i trrcit to ILAST, 

^cS^ D (U,S r F. I ^S n SS t f S- ™f the moduius 
caU a SS y ourve A fit?in| “broutinfLsSF (passing XFREQ) 

to'calculate'the intelpolated f regency index FO P^and 

the modulus ?^n call'^F again (pass- 

1P dlflf LSCF°a f tSrd“?me' (passing^SHALfS tinT**' 

l h ? R IS a i ^ FQ 5S?“if r^ula^rihfn 
1S th index of the frequency search band mid- 

point, XMAX = CABS (AWO (KF) ) , PHASEI = AIMAGI ( AWO ( (KF .) ) , 

SI 

period, subtract one period f ro “ IT ^' J'-* ITMAX - 
bT is greater than (1.0/FREQ), then ITMAX - ITMAX 
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I "T( 1 i O/(FREQ^DT ) ). Print to^ejile^ndicated 

linen? unle?? 1 FREQ* i^nter^than^O.0035, f in^hich 

EkSSs: sks2 

and AIMAG ( AWO ( I ) ) to the output file. 


Subroutine HANN (LA, All, BIAS) 


Definition of Variables 

Variables passed as arguments (equals NTB1) 

“ ‘-"SST. “iS -nn‘w?ndow is ap- 

BIAS - P mean of array All (after windowing) 

Constants 

PI — set to 3.1415926 

^-"aSa^hoMinc, the calculated discrete Hann window 

ITM - index of midpoint of array All 
RM - ITM as a real variable 

Subroutine Calls 


MEAN (LA, All, BIAS) 
Code Discussion 


The array All, and LA, the length of All, ^re 

m°: «2". tonve ted^o n r i eal a Value d PH. iF or 

ffy? IX SScil ™.t * (iL - COS (PI* ITI/RM), 

and a PP ly J*® i Han 2ali n sSbr^tine MEAl^to find' the mean 
value } BIA^ of the C windowed°array All, and subtract the _ 
windowed BIAS from the windowed array, All(I) All( ) 
HW (I) * BIAS. 
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Subroutine MEAN (LA, A22, SA) 


Definition of Variables 


Variables passed as arguments 
LA - length of the data array A22 
A22 - data array 
SA - mean value of array A22 

Code Discussion 


The values of the array A22, and LA, the length of 
array A22, are passed by the calling statement. SA is 
the sum of the values of the individual entries of A22, 
divided by LA. 


Subroutine F0UR1 (DATA, NN, ISIGN) 

Subroutine F0UR1 is the standard FFT routine found 
in "Numerical Recipes". 

Definition of Variables 


Variables passed as arguments 

DATA - data array (complex, but converted to a real 
array of double length) 

NN - number of points in the Fourier Transform 
ISIGN - +1 indicates forward transform, -1 inverse 

Other Variables 

WR, WI, WPR, WPI, WTEMP, THETA - all double precision 
variables used in the usual array shuffling procedures 

Code Discussion 


Subroutine F0UR1 utilizes the array shuffling pro- 
cedure common to most FFT routines. For more details, 
consult "Numerical Recipes", or other sources that dis- 
cuss FFT routines at length. 
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Subroutine LSCF (PO, FMAX, U_IN, IOPT, G_FLAG ) 


Definition of Variables/Flags 


Variables passed as arguments 
PO - interpolated location of maximum 
FMAX - value of the maximum located at PO 
U_IN - input ordinate array (length of 7) 

IOPT - action option, either 1 or 2 

IOPT = 1: Find PO and compute maximum at PO 
IOPT = 2 : Only compute value at PO 
G_FLAG - indicates data validity, TRUE if the peak is 
inside the search zone, FALSE if outside 

Other Variables 

US17 - sum of U_IN ( 1) and U_IN(7) 

US3 5 - sum of U_IN(3) and U_IN(5) 

COF1 - constant term in fitting quadratic 

COF2 - coefficient of linear term in fitting quadratic 

C0F3 - coefficient of square term in fitting quadratic 

Code Discussion 


Subroutine LSCF does a least squares curve fit of 
a quadratic to 7 data points sampled at integral inter- 
vals. The indices of the 7 points range from -3 to +3. 
The quadratic is F(P) = C0F1 + COF2*P + COF3*P**2, the 
max occurs at PO = -C0F2/ (2*COF3) , and the maximum 
value is F(PO) = COF1 - (COF2**2) / (4 . 0*COF3 ) . The 3 
coefficients are computed by multiplying the 3x7 matrix 
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times the array U_IN (7 entries in the array) , with the 
results C0F1 = first row x U_IN, COF2 = second row x 
U_IN, and C0F3 = third row x U_IN. The array U^IN and 
the option variable IOPT are passed by the callTng 
statement. To further use the symmetry of the matrix, 
US17 = U_IN ( 1 ) + U_IN ( 7 ) and US35 = U_IN(3) + U_IN(5) 
are created. C0F1, COF2 , and COF3 are calculated as 
described above. C0F3 is checked to be sure that it 
does not equal 0.0 (equaling 0.0 would prevent calcula- 
tion of the maximum F(P0) = COF1- (COF2**2) / (4.0*COF3) ) ; 
if COF3 = 0.0 then G_FLAG is set to false and control 
returns to the calling subroutine WORK. 

For IOPT option 1, calculate P0 = -0. 5*COF2/COF3 , 
FMAX = (C0F1 + 0.5*P0*COF2) /84, and set G_FLAG to true. 
To check if the calculated P0 is valid, compare the 
absolute value of P0 to 3; if ABS(PO) greater than 3, 
then set G_FLAG to false. For IOPT option 2, calculate 
FMAX = (C0F1 + P0* (C0F2 + PO*COF3 ) ) /84 . 
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Subroutine FBAND (FL, FH, TLKM, PF) 


Definition of Variables 


Variables passed as arguments 

FL - low frequency bound of the frequency search band 
PH _ high frequency bound of the frequency search band 
TLKM - tether length in kilometers 

pp — percentage used to compute frequency search band 
Common Variables 

DENSITY - tether density in kg/km (set to 8.35 kg/km) 
TOTALL - total tether length (set to 22.0 km) 

MSAT - satellite mass in kg (set to 510.0 kg) 

MORB - orbiter mass in kg (set to 100000.0 kg) 

ALTKM - orbit altitude in km (set to 325.0 km) 

Other Variables _ 

FC - center frequency of the frequency search band - 
0.5 * SQRT ( CK * MSTAR/TLKM) 

(estimate of the skiprope frequency based on the 
tether parameters listed as common variables and 
the average tether length TLKM) 

OMSQ - orbit rate squared (OMSQ = ORBRAT ESQ (ALTKM) ) 

CK - working variable = 3.0 * OMSQ / DENSITY 
MO - sum of orbiter mass and tether mass (but not the 
satellite mass!) MO = MORB + TOTALL * DENSITY 
Q - working variable =0.5 * DENSITY * TLKM 
MSTAR - working variable = ( (MO-Q) * (MSAT+Q) ) / (MO+MSAT) 
DP - fraction of FC to be subtracted from FC to create 
FL and added to FC to create FH (DF = FC*PF/100.0) 

Function Call 


REAL FUNCTION ORBRATESQ (ALTKM) 

Parameters . , . , 

GM - acceleration due to gravity, in meters/sec**2 (set 

to 9.81098) 

PE - radius of earth in km (set to 6378.17) 


Other Variables 
R - GM/ (1000.0 * RE) 

ORBRATESQ - R/(1.0 + ALTKM/RE) **3 

Code Discussion 


The values of the average tether length in km, 
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TLKM, and the percentage around the estimated skiprope 
frequency, PF, are passed as arguments of the calling 
statement, and tether/orbiter parameters DENSITY, 

TOTALL, MSAT, MORB, and ALTKM are passed as common 
variables. Using equations derived from the dynamical 
analysis of the skiprope frequency vs tether length, 
the estimated skiprope frequency FC is calculated (see 
the variable list for the equations used) . Search band 
bounds FL and FH are computed from FC by using the per- 
centage PF, FL = FC-FC*PF/100 . 0 and FH = FC+FC*PF/100 . 0 . 


Subroutine READINDATA (TO ,TF, X, Y,TLG,MODE, LF ,MFLAG) 


Definition of Variables/Flags 


Variables passed as arguments 
TO - time tag for first point in buffer 
TF - time tag for last point in buffer 
X(I) - x axis gyro data array (deg/sec) 

Y ( I ) - y axis gyro data array (deg/sec) 

TLG(I) - tether length array (kilometers) 

MODE - integer variable, signifying the amcsmode for 
the last time point 'TF' 
amcsmode = 0 indicates no valid data 
amcsmode = 1 indicates passive case 
amcsmode = 2 indicates yaw hold 
amcsmode = 3 indicates spin case 
LF - number of time points stored in each array 
MFLAG - logical flag to denote if the amcsmode 

changed from 1 or 2 to a 0 or 3 between the 
times of TO and TF 

Other Variables 

JMODE - amcsmode of time tag TO 
External Input File 


IDFTXY.DAT 

The external input file IDFTXY.DAT simulates the 
preprocessed telemetry data stream. Each line has the 
following format: 

TIME, X(I), Y ( I ) , TLG(I), MODE 

Code Discussion 


The external input file IDFTXY.DAT is opened as 
logical unit 10 and the first line read, with the time 
recorded as TO and the MODE value as JMODE. A read 
loop is entered, and lines will be read as long as the 
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MODE remains unchanged or if the MODE only changes from 
1 to 2 or 2 to 1, until the end of the file. If the 
MODE changes other than from 1 to 2 or 2 to 1, then 
reading stops, and MFLAG is set to false. The last 
line read, whether or not the end of file is reached, 
has the time recorded as TF, the MODE value returned to 
the main program, and is the point where LF is calcula- 
ted. The file IDFTXY.DAT is then closed. 


Subroutine ODTF (TO, TM, TF, DT, LE, LB, LEB, TL, AX, 

FX, PX, AY, FY, PY, FA, TA, WK, PSIGN, 
UMAX, VMAX, FL, FH, AVGX, AVGY) 


Definition of Variables 


Variables passed as arguments 

TO - first time point in buffer (mission elapsed time) 
TM -time at midpoint of data window (met) 

TF - last time point in buffer (met) 

DT - sample time (sec) 

LE - index of last point in data window 

LB - index of first point in data window 

LEB - number of points used in data window 

TL - tether length used for this data window (km) 

AX - peak magnitude of x axis gyro rate (deg/sec) 

FX - calculated skiprope frequency in x axis (hz) 

PX - phase angle in x axis 

AY - peak magnitude of y axis gyro rate (deg/sec) 

FY - calculated skiprope frequency in y axis (hz) 

PY - phase angle in y axis 

FA - average frequency with G_FLAG set to true 
TA - average period, reciprocal of FA 

WK - conversion factor from rate data to skiprope amli- 
tude (meter-sec/deg) 

PSIGN - skiprope polarity w.r.t Z LVLH axis 
UMAX - maximum in plane midnode skiprope amplitude 
VMAX - maximum out of plane midnode skiprope amplitude 
FL - lower boundary of the frequency search band (hz) 

FH - upper boundary of the frequency search band (hz) 
AVGX - computed mean of the Hann-windowed x axis gyro 
rate data 

AVGY - computed mean of the Hann-windowed y axis gyro 
rate data 

Output File 


F_RECORD.DAT (with logical unit number 18) 

The output file F_RECORD.DAT has 6 lines with the 
following format: 
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TM 

LB, 

LE, 

LEB 



TO, 

TF, 

DT, 

TL, WK, 

PSIGN 

AX, 

FX, 

PX, 

AY, FY, 

PY 

FA, TA, UMAX, VMAX, 
AVGX, AVGY 

FL, FH 


Code Discussion 


The variables in the variable list are passed as 
arguments and written to logical unit 18 (output file 
F RECORD.DAT) in the format described above. 
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Quick Reference 


Required External Files 


PARAM.DAT (read in the main program) 

The external input file PARAM.DAT consists of two lines 
having the following format: 

PF SRCH_BAND, R_ARM, ODF_TIME, ODF_FFT 
DENSITY, TOTALL, MS AT, MORB, ALTKM 

IDFTXY.DAT (read in the subroutine READINDATA) 

The external input file IDFTXY.DAT simulates the pre- 
processed telemetry data stream. It consists of a 
maximum of 3000 lines with the following format: 

TIME, X(I), Y ( I ) , TLG(I), MODE 


Operator Input Prompts (in order of appearance in the 
‘ main program) 

Prompt: Yaw maneuver calculations will not be performed 
unless requested. Type yes (y) to calculate, 
no (n) otherwise. 

Prompt: DT is (DT value) , is this okay to use - yes (y) 
or no (n) 

If no, then prompt: Enter DT in seconds 

Prompt: Recommend using at least (LEAST value) points 
for 3 data cycles. Last point in buffer is 
(LF) . What is the starting time index? 

After receiving the start time index LB, the prompt 
continues with: What is the last time index? 

Prompt: Tether length is ( TLNGTH value) kilometers - ok 
to use, reply with a yes (y) or a no (n) 

If no, then prompt: Enter the tether length in kilome- 
ters, and then repeat original tether length prompt. 

Prompt: FLOW & FHIGH = (FLOW, FHIGH values) 

Are these bounds okay to use - y or n 
If no, then prompt: Enter the two values in hz, and 
then repeat original frequency bounds prompt. 

Prompt: Enter orbiter yaw axis angle in degrees 

Verification prompt: Orbiter nose is degrees wrt 

X-LVLH axis. Is this correct (Y or N) ? 

If no, repeat the original orbiter yaw axis angle 

prompt. 
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Warnings/Abort Situations 


If mode flag M_FLAG is false, the following warning is 
printed to the screen: 
************************************* 

WARNING 

************************************* 

AMCSMODE STATUS CHANGED DURING THIS DATA STREAM 
NUMBER OF POINTS REDUCED - DATA SET HAS SAME MODE 
WILL SET AMCSMODE TO LAST READING 

If AMCSMODE = 0, then print this warning: 
************************************************* 
AMCSMODE INDICATES NO VALID DATA - PROGRAM ABORTS 
************************************************* 

The following warning is printed if AMCSMODE = 3 
********************************************* 

AMCSMODE INDICATES SPIN CASE - PROGRAM ABORTS 
********************************************* 

Warning for G_FLAGY false and G_FLAGX true: 

Y AXIS DATA IS SUSPECT - FREQ OUT OF BAND 

FREQUENCY RETURNED IS PREDICTED MIDPOINT 

DATA SHOULD NOT BE USED WITHOUT CAUTION 

NEITHER POLARITY NOR YAW MANEUVER CALCULATIONS ARE 

PERFORMED 

Warning for G_FLAGX false and G_FLAGY true: 

X AXIS DATA IS SUSPECT - FREQ OUT OF BAND 

FREQUENCY RETURNED IS PREDICTED MIDPOINT 

DATA SHOULD NOT BE USED WITHOUT CAUTION 

NEITHER POLARITY NOR YAW MANEUVER CALCULATIONS ARE 

PERFORMED 

Warning for both G_FLAGX and G_FLAGY false: 

BOTH AXES ARE BAD ***** 3 SUGGESTIONS 

1) SUGGEST LOOK AT TIME PLOTS OF GYRO SIGNALS. 

2) SUGGEST MAKE FFT PLOTS AND LOOK AT DATA. 

3) SUGGEST WIDENING SEARCH BAND. 

If the curve fitting subroutine LSCF calculates the 
quadratic coefficient COF3 = 0.0, print this warning: 
****************** WARNING ******************* 
QUADRATIC COEFFICIENT EQUALS ZERO - CANNOT COMPUTE A 
MAXIMUM FREQUENCY VALUE. 
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INTRODUCTION 


This document describes the Test Plan and Procedures for evaluating the 
Frequency Domain Skiprope Observer. The plan is divided into two parts, one 
for Station 2 conditions, and the second for Station 1 conditions. No concrete 
performance requirements exist at Station 1, because the main focus of the 
Observer is to support a yaw maneuver at Station 2. Nevertheless, two tests 
using simulations at Station 1 are included to demonstrate Observer 
performance. Additional Station 1 tests using model test signals are also 
included to demonstrate and/or define Observer performance boundaries at 
Station 1. 

The test cases enumerated in the ECR, using both model gyro signal plus 
actual simulation data, will be fully documented with input data and filter 
output results. These cases should be used to verify observer code whenever 
the code is transferred to a different computer system. 

This Test Plan calls for using simulated gyro noise. The noise source to 
be used for all testing is a portable random number generator as documented in 

Reference and included as Appendix 2. A herein. The model for 

generating a Gaussian distribution for noise is also a part of Appendix 2. A. 


PART I - STATION 2 CONDITIONS: (2.4 km TETHER LENGTH) 


This part is divided into four Test Groups: 

A. Six cases using model gyro signals per ECR 

B. Three cases using simulation data per ECR 

C. Systematic error testing without noise (using model gyro signals) 

D. Systematic error testing with noise (using model gyro signals) 


The purpose of the first two groups is to establish test case results for 
code transfer as well as prove performance of the Observer. Where noise is 
modelled using the random number generator, documentation of these cases will 
include the initial value of the random number seed. Users of the test plan 
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-should use any or all of the fully documented cases to verify the code - 
results should vary from the results documented only to within expected 
roundoff error. Users should also vary the initial value of the random number 
seed to fully statistically test the Observer. 

Cases in Group A will verify the ECR requirements in paragraph 3.2. (1) 
and 3 . 2 . (2), i.e., the error of the angular rate amplitude shall not exceed 2%, 

“with a maximum total phase error of 25 degrees after 15 minutes from the last 
time point used in the data window. 

_ Cases 1 and 2 in Group B use simulation data to verify that the secondary 

Observer meets the performance requirements for the primary Observer. 

Amplitude and phase errors will be measured on a root-mean-square basis as 
outlined in the ECR, paragraphs 5.3.2 and 5.3.1. 

Case 3 of Group B will be tested the same as cases 1 and 2; however, the 
conditions inherent to this case violate the constraints and limitations 

- enumerated in the ECR. This case is included to demonstrate trends in 
degradation for highly transient conditions. 

__ Cases in Group C are used to best define the ultimate performance of the 

Observer under ideal (noise-free) conditions. Model gyro signals are inputted 
to the Observer for a range of frequencies at Station 2 and for various 
pairings of skiprope and pendulous phases. Percentage errors between the input 

“values and output values of the skiprope amplitude, frequency, and phase are 
reported. 

_ Cases in Group D are used to define the expected performance in a noisy 

environment. For each case in Group C, 50 different noisy signals are 
generated (using the portable random number generator and the Box-Muller 
algorithm to generate Gaussian distributed noise described in Appendix 2. A). 

“The maximum amplitude, frequency, and phase errors found in the set of 50 noise 
runs are reported, as well as the average errors over the 50 runs. For each 
frequency tested (11 total), the largest maximum and largest average errors are 

— also reported. 
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A: STANDARDIZED TEST CASES 


These test cases are specified in the ECR and have a model gyro signal of the 
following form: 

N (T) + AO + Al*COS(2*PI*Fl*(I-l) *DT + PHI1) + A2*COS (2*PI*F2 * ( 1-1) *DT + PHI2) 

where N(T) = GAUSSIAN DISTRIBUTED NOISE WITH SIGMA = 2.8E-04 DEG/S 
AO = ORB RATE (DEG/S) 

A1 = SKIPROPE AMPLITUDE (DEG/S) 

FI = SKIPROPE FREQUENCY (HZ) 

PHIl = SKIPROPE PHASE (DEG) 

A2 = PENDULOUS AMPLITUDE (DEG/S) 

F2 = PENDULOUS FREQUENCY (HZ) 

PHI2 = PENDULOUS PHASE (DEG) 

DT = SAMPLE TIME (SEC) 

Input data files to the Observer must be in the following format: 

TIME, X GYRO SIGNAL, Y GYRO SIGNAL, TETHER LENGTH, AMCSMODE 

Both the x gyro signal and y gyro signal are of the form detailed above. 

Please note that the x and y gyro signals should be generated independently, 
albeit concurrently, and have distinct Gaussian distributed noise terms (as 
detailed explicitly in the discussion of the noise generation found in Appendix 
2. A). At Station 2 the tether length is 2.4 km, and amcsmode can be either 1 
or 2 ! The time values must be spaced at the sample rate, preferably 1.024 s. 
Users may write their own signal generating programs, or they may use the 
program CREATE. FOR listed in Appendix 2. A. 

The following table is the list of values used to generate the model gyro 
signals. Please note the following: 

1) The skiprope amplitude and frequency values, and the pendulous amplitude 
values are used for both the x and y signals. 

2) The orb rate is added only to the y gyro signal (AO = 0.0 for the x gyro 
signal) . 

2 j Skiprope phase is the y axis value. The x axis value is phase + 9 0.0. 

4) Pendulous phase is the x axis value. The y axis value is phase - 90.0. 

5) Gaussian distributed standard deviation noise of 2.8E— 03 deg/s is 10 times 
the expected noise sigma of 2.8E-04 deg/s. 

6) All test cases should use data lengths of at least 3 skiprope periods and a 
sample time of 1.000 s or 1.024 s (one skiprope period is calculated as 
1.0 / (freq x sample time) ). 

7) Pendulous frequency = 0.03125 Hz 
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8) Cases with non-zero noise should use 10 runs each. The averages of the 10 
runs constitute the output. 


CASE 

NO. 

NOISE 

ORB 

RATE 

AMP 

SKIPROPE 

FREQ 

PHASE 

PENDULOUS 
AMP PHASE 

1 

0.0 

0.06 

0.02 

0.0054 

163.0 

0.0 

10.0 

2 

2.8E-03 

0.06 

0.02 

0.0054 

163.0 

0.0 

10.0 

3 

0.0 

0.06 

0.02 

0.0046 

-40.0 

0.5 

60.0 

4 

2.8E-03 

0.06 

0.02 

0.0046 

-40.0 

0.5 

60.0 

5 

2.8E-03 

0.06 

0.15 

0.0054 

70.0 

0.0 

50.0 

6 

2.8E-03 

0.06 

0.15 

0.0054 

70.0 

0.5 

50.0 


Appendix 2.B lists tables of results for all six cases and 40 initial values of 
the random number seed (4 test cases with non-zero noise x 10 noise runs each) . 


B: SIMULATION RUNS (VERIFICATION MATRIX) 


For this group of the simulations, the user should use the option to print the 
calculated rates in the Observer program (UNOMSC.FOR) ,i.e. , set ODF_TIME to 
true. A root_mean_square comparison of the skiprope amplitudes and phases is 
then performed between the original simulation data and the data generated by 
the Observer. A sample comparison program is listed in Appendix 2.C. 


CASE SKIPROPE (IN PLANE X OUT OF PLANE) TETHER LENGTH (km) 


1 

20 

X 

20 

2.4 

2 

60 

X 

60 

2.4 

3 

80 

X 

40 

2.4 


Notes: 

1) Data length should be at least 3 skiprope periods with a sample time of 
1.000 s or 1.024 s. 

2) These 3 simulations are required by the ECR. The observer should work 
properly given any valid simulation, i.e., a simulation without spin. 
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c- SYSTEMATIC ERROR TESTING ACROSS A RANGE OF SKIPROPE 
FREQUENCIES AND SKIPROPE / PENDULOUS PHASE PAIRINGS 
NOISE-FREE CASE 


Model qvro signals are of the same form as detailed in section A. These tests 
a?e of the gylo signal itself, so only one axis is necessary (all parameters 
-are used on this one axis). The programs listed in Appendix 2. D automate the 
procedure of creating the signals and running the essentials of the Observer by 
incorporating various loops into the body of the program to eliminate user 
_ input in creating data files and/or running the Observer. 


For each frequency in the range (0.0045 - 0.0055 Hz), at intervals of 0.0001 Hz 
(a total of 11 frequencies), run the following tests: 


Parameters : 


_ orb rate = 0.065 deg/s 
' skiprope amp = 0.02 deg/s 

pendulous frequency = 0.03125 Hz 
pendulous amplitude = 0.5 deg/s 
“data length = at least three periods of data 
noise = 2.8E-03 deg/s 
sample rate = 1.024 s or 1.000 s 


1} Varv the skiprope phase from -180.0 deg to 180.0 deg in increments of 10 
deq For each skiprope phase vary the pendulous phase from -180.0 deg to 
180 ! 0 deg in increments of 10 deg (a total of 37 x 37 = 1369 cases). 


2) In each case record the per cent errors in the calculated skiprope 
amplitude, frequency, and phase. 

3) Find the maximum amplitude, frequency, and phase per cent errors and the 
associated skiprope and pendulous phases for each. 


4) Plot error surfaces of the amplitude, frequency, and phase errors (a total 
of 33 plots - 3 plots for each frequency x 11 frequencies) . 


D* SYSTEMATIC ERROR TESTING ACROSS A RANGE OF SKIPROPE 
FREQUENCIES AND SKIPROPE / PENDULOUS PHASE PAIRINGS 
NOISE CASE 


Model gyro signals are of the same form as detailed in section A. These tests 
are of 9 the gyro signal itself, so only one axis is necessary (all parameters 
are used on this one axis). The programs listed in Appendix 2. D automate the 
procedure of creating the signals and running the essentials of the Observer by 
incorporating various loops into the body of the program to eliminate user 
input in creating data files and/or running the Observer. 
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For each frequency in the range (0.0045 - 0.0055 Hz), at intervals of 0.0001 Hz 
~ (a total of 11 frequencies) , run the following tests: 

Parameters : 


orb rate = 0.065 deg/s 
skiprope amp = 0.02 deg/s 
—pendulous frequency = 0.03125 Hz 
pendulous amplitude = 0.5 deg/s 
data length = at least three periods of data 
noise = 2.8E-03 deg/s 
— sample rate = 1.024 s or 1.000 s 

1) For each of the 1369 phase relationship cases listed in C.l for the 
noise-free case, run 50 noise runs (Gaussian distributed noise) . 

2) In each case record the average and maximum per cent errors in the 
_ calculated skiprope amplitude, frequency, and phase. 

3) Find the maximum per cent errors in the amplitude, frequency, and phase and 
the associated skiprope and pendulous phases for each. 

4) Find the largest maximum errors in the amplitude, frequency, and phase and 
the associated skiprope and pendulous phases for each. 

5) Plot representative samples of the surfaces generated m part b) . 


PART II - STATION 1 CONDITIONS: (20.0 km TETHER LENGTH) 


This part is divided into three Test Groups: 


A. Two cases using simulation data per ECR 

B. Systematic error testing without libration component — . both noise— free and 
noise cases (using model gyro signals) 


C. Systematic error testing with libration component - both noise-free and 
noise cases (using model gyro signals) 

Cases 4 and 5 in Group A use simulation data to verify that the secondary 
Observer meets the performance requirements for the primary Observer. 

Amplitude and phase errors will be measured on a root-mean-square basis as 
outlined in the ECR, paragraphs 5.3.2 and 5.3.1. (Note: The numbering of the 
simulation cases follows the convention of the ECR, which has simulations from 
both stations 2 and 1 in one table and numbered sequentially - cases 1, 2, and 
3 at station 2 cases 4 and 5 at station 1.) 
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Cases in Group B are used to define the expected performance of the 
Observer at Station 1 without the influence of a libration component. Both the 
noisy and ideal (noise-free) environments are considered. 

Cases in Group C are used to define the expected performance of the 
Observer at Station 1 with the influence of a libration component. Both the 
noisy and ideal (noise-free) environments are considered. 

For both Groups B and C, runs are performed for all cases with data 
lengths of 2, 3, and 4 skiprope periods. Noise runs are performed with both 
1 sigma and 3 sigma Gaussian distributed noise standard deviations. 


A: SIMULATION RUNS (VERIFICATION MATRIX) 


For this qroup of the simulations, the user should use the option to print the 
calculated rates in the Observer program (UNOMSC.FOR) ,i.e. , set ODF_TIME to 
true. A root mean square comparison of the skiprope amplitudes and phases is 
then performed between the original simulation data and the data generated by 
the Observer. A sample comparison program is listed in Appendix 2.C. 


CASE 


SKIPROPE (IN PLANE X OUT OF PLANE) 


TETHER LENGTH (km) 


4 

5 


80 

80 


40 

80 


20.0 

20.0 


No^cs • • 

1) Data length should be at least 3 skiprope periods with a sample time of 

1.000 s or 1.024 s. , _ . ,, . . 

2) These 2 tests are required by the ECR. The observer should be able to work 
properly given any valid simulation, i.e., a simulation without spin. 


B: SYSTEMATIC ERROR TESTING OF MODEL SIGNALS WITHOUT LIBRATION 
USING DATA LENGTHS OF 2 , 3, OR 4 SKIPROPE PERIODS 


Model gyro signals are of the same form as detailed in section A, part I. The 
CREATE FOR program listed in Appendix 2. A can be used to generate signals or 
this group. Appendix 2.D lists programs that automate the data file generation 
and Observer testing for the cases in this group. 
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?or each of the desired data lengths, data files should be generated with the 
following values: 


PARAMETERS 

• 

« 


LIBRATION 

SKIPROPE 

PENDULOUS 

* 

AMPLITUDE 

(X 

AXIS) 

0.004 

0.0034 

0.05 

(units = deg/s) 

— 

(Y 

AXIS) 

0.004 

0.0034 

0.05 

FREQUENCY 

(X 

AXIS) 

1/2713 

0.0019 

0.089 

(units = Hz) 


(Y 

AXIS) 

1/3132 

0.0019 

0.089 

PHASE 

(X 

AXIS) 

varies 

varies 

0.0 

(units = deg) 

SIGMA = 2. 

(Y 

8E- 

AXIS) 

-04 deg/s 

varies 

varies 

-90.0 


3RB RATE = 0.065 deg/s 
_DT = 1.024 s or 1.000 s 

1) For the model signal without libration component vary the skiprope phase 
from -180.0 deg to 180.0 deg in increments of 10 deg (a total of 37 cases). 

2) In each case record the per cent errors in the calculated skiprope 
amplitude, frequency, and phase. 

3) Find the maximum amplitude, frequency, and phase per cent errors and the 
associated skiprope phases for each. 

~4) Plot error curves of the amplitude, frequency, and phase errors (a total of 
3 plots) . 

-Oo the following steps for Gaussian distributed noise signals, using 

noise = 1 x sigma = 2.8E-04 and noise = 3 x sigma = 8.4E-03: 

5 ) For each of the 37 phase relationship cases listed in 1) for the noise-free 
case, run 50 noise runs (Gaussian distributed noise) . 

5) In each case record the average and maximum per cent errors in the 
— calculated skiprope amplitude, frequency, and phase. 

7) Plot representative samples of the curves generated in part 6) . 


C: SYSTEMATIC ERROR TESTING OF MODEL SIGNALS WITH LIBRATION 
USING DATA LENGTHS OF 2 , 3, OR 4 SKIPROPE PERIODS 


Model gyro signals are of the same form as detailed in section A, part I, with 
the addition of a libration component term of the form: 


ALIB*COS (2*PI*FLIB* (1-1) *DT + PHLIB) 


ALIB = LIBRATION AMPLITUDE 
FLIB = LIBRATION FREQUENCY 
DT = SAMPLE RATE 
PHLIB = LIBRATION PHASE 


ADDendix 2. A lists a program (CRELIBR. FOR) that generates a model gyro signal 
with a libration component. Appendix 2.D lists programs that automate the 

generation and Observer evaluation for the cases in this group. 


For each of the desired data lengths, data files should be generated with the 
following values: 


PARAMETERS : 


AMPLITUDE (X AXIS) 
(Y AXIS) 

FREQUENCY (X AXIS) 
(Y AXIS) 

PHASE (X AXIS) 

(Y AXIS) 


LIBRATION 

SKIPROPE 

PENDULOUS 


0.004 

0.004 

0.0034 

0.0034 

0.05 

0.05 

(units. = deg/s) 

1/2713 

1/3132 

0.0019 

0.0019 

0.089 

0.089 

(units = Hz) 

varies 

varies 

varies 

varies 

0.0 

-90.0 

(units = deg) 


SIGMA = 2.8E-04 deg/s 
ORB RATE = 0.065 deg/s 
DT = 1.024 s or 1.000 s 


H For the model signal with libration component vary the skiprope phase from 
-180.0 deq to 180.0 deg in increments of 10 deg. For each skiprope phase, 
vary the libration phase from -180.0 deg to 180.0 deg (a total of 1369 cases 
37 skiprope phases x 37 libration phases) . 

2) In each case record the per cent errors in the calculated skiprope 
amplitude, frequency, and phase. 

3) Find the maximum amplitude, frequency, and phase per cent errors and the 
associated skiprope and libration phases for each. 

4) Plot error surfaces of the amplitude, frequency, and phase errors (a total 
of 3 plots) . 
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Do the following steps for Gaussian distributed noise signals, using 

noise = 1 x sigma = 2.8E-04 and noise = 3 x sigma = 8.4E-03: 

5) For each of the 1369 phase relationship cases listed in 1) run 50 noise runs 
(Gaussian distributed noise) . 

6) In each case record the average and maximum per cent errors in the 
calculated skiprope amplitude, frequency, and phase. 

7) Find the maximum average per cent errors in the amplitude, frequency, and 
phase and the associated skiprope and libration phases for each. 

8) Find the largest maximum errors in the amplitude, frequency, and phase and 
the associated skiprope and libration phases for each. 

9) Plot representative samples of the surfaces generated in part 6) . 
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APPENDIX 1 



noon onooo 


c 

c 

c 

c 


NAME IS UNOMSC . FOR (BACK UP SKIPROPE OBSERVER) 
THIS VERSION IS COMBINED FROM UNO AND MSFC 
DATE IS SEPTEMBER 1991 


INTEGER NDIM, NFT, NPNT, AMCSMODE 
PARAMETER (NDIM= 3000, NFT=8192, NPNT=7) 

PARAMETER (PI=3 . 1415926 , DFR=57 . 2957795 , RFD-0 . 017453292 ) 
REAL* 4 X (NDIM) , Y (NDIM) , T LG (NDIM) 

CHARACTER* 1 REPLY 
CHARACTER*8 POLARITY 

logical g_flagx, g_flagy, s_flag, m_flag, odf_time, 
0DF-FF^ ealA4 pFJ3RCHJ3AND> r _arm, density, totall, msat, morb, 

ALTKM C0MM0N ^ FREQ ^ DENSITYf TOTALL, MSAT, MORB, ALTKM 

COMMON lb, le, dt, df, flow, fhigh 

C READ IN DATABASE PARAMETERS FROM FILE ' PARAM . DAT ' 

C PF SRCH BAND IS % NUMBER TO COMPUTE SEARCH BAND. 

R ARM IS DISTANCE FROM ORBITER C.M. TO CENTER LINE 

” OF DEPLOYER BOOM. „ m __ 

ODF TIME IS LOGICAL FLAG TO CONTROL PRINTING OF TIME 
DATA TO OUTPUT FILE. THIS FLAG IS NOMINAL .FALSE. 
MEANING TIME DATA IS NOT PRINTED. 

ODF FFT IS LOGICAL FLAG TO CONTROL PRINTING OF FREQUENCY 
DOMAIN DATA (FFT'S) TO OUTPUT FILE. THIS FLAG IS 
NOMINAL .TRUE. .. THE DATA IS PRINTED OUT. 

DENSITY IS TETHER DENSITY IN KG PER KM — 8.35 KG/KM. 
TOTALL IS TOTAL TETHER LENGTH = 22.0 KILOMETERS. 

MSAT IS SATELLITE MASS IN KGS. DEFAULT = 510. 

MORB IS ORBITER MASS IN KGS. DEFAULT = 100,000. 

ALTKM IS ORBIT ALTITUDE IN KM. DEFAULT = 325. KM. 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OPEN ( 10 , FILE= ' PARAM . DAT ' , STATUS= ' OLD ' ) 

READ (10,*) PF SRCH BAND, R_ARM, ODFJTIME, ODF_FFT 
READ (10,*) DENSITY, TOTALL, MSAT, MORB, ALTKM 
CLOSE ( 10 , STATUS= ' KEEP ' ) 

OPEN FILE COMMANDS 

THESE ARE OUTPUT FILES FOR RECORD. 

CLOSE STATEMENTS APPEAR JUST PRIOR TO 'END' STATEMENT. 

OPEN ( 11 , FILE=' F FFTX.DAT' ,STATUS=' UNKNOWN' ) 

OPEN ( 12 , FILE= ' F FFT Y . DAT ' , STATUS=' UNKNOWN' ) 

OPEN ( 13 , FILE= ' F TXYUV.DAT' , STATUS=' UNKNOWN ' ) 

OPEN ( 17 FILE—' F YAWMAN . DAT' , STATUS= ' UNKNOWN ' .) 

OPEN ( 18 , FILE= ' F~RECORD . DAT ' , STATUS= ' UNKNOWN ' ) 

THE RATE GYRO DATA SHOULD ALWAYS BE IN LVLH FRAME AND 
IS THE RATE RELATIVE TO LVLH. I.E. ORBITAL RATE HAS 
ALREADY BEEN REMOVED. NOTE: PROGRAM WILL STILL WORK 
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c 

c 

c 

c 


PROPERLY IF ORBITAL RATE IS NOT REMOVED. 

GO READ FILE FOR TELEMETRY DATA 

CALL READINDATA (TO, TF , X, Y, TLG, JMODE, LF, M_FLAG) 

CHECK ON M FLAG STATUS . . . SET AMCSMODE BY M_FLAG AND JMODE . 


PRINT*, 'ESTIMATED TETHER LENGTH IN KILOMETERS IS ' , TLG (LF) 
IF (M_FLAG) THEN 
AMCSMODE = JMODE 

C AMCSMODE STATUS ON TELEMETRY NEVER CHANGED FOR ALL 

POINTS, 

C 

ELSE 

PRINT * ' *************************************' 

PRINT * , ' WARNING 1 

PRINT * , ' *************************************' 

PRINT *,' AMCSMODE STATUS CHANGED DURING THIS DATA 


OR CHANGED BETWEEN MODES 1 & 2 ONLY. 


STREAM' 
MODE . ' 


PRINT * , ' NUMBER OF POINTS REDUCED - DATA SET HAS SAME 


PRINT *,' WILL SET AMCSMODE TO LAST READING' 

AMCSMODE = JMODE 
END IF 

PRINT *,' AMCSMODE IS = ', AMCSMODE 
IF (AMCSMODE. EQ. 0. ) THEN 
PRINT* 

PRINT* f ************************************************* ' 

' PRINT*, 'AMCSMODE INDICATES NO VALID DATA - PROGRAM 

ABORTS ' 

PRINT *, ’*************************************************' 

CALL EXIT 
ELSE 

IF (AMCSMODE. EQ. 3) THEN 
PRINT* 

PRINT* / ********************************************* , 

' PRINT*, 'AMCSMODE INDICATES SPIN CASE - PROGRAM 

ABORTS ' 

PRINT * , ********************************************** ’ 

CALL EXIT 
END IF 
END IF 

C SET VALUE FOR FLAG (S FLAG) TO CONTROL COMPUTATIONS OF YAW 

C MANEUVER OUTPUTS. THIS IS BASED ON AMCSMODE AND TETHER 

C E NGTHq PERAT0R FULL CONTROL OF OPTION TO EXECUTE THIS 

COMPUT AT I pLAG (logical) USED TO DENOTE IF YAW MANEUVER 
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COMPUTATIONS 

C ARE REQUIRED. IF ' TRUE ' MEANS DO THE YAW 

COMPUTATIONS. 

c IF 'FALSE' MEANS DO NOT DO THE 

COMPUTATIONS. 

C 

S FLAG = .FALSE. 

PRINT* , 'YAW MANEUVER CALCULATIONS WILL NOT BE PERFORMED 
UNLESS ' 

PRINT*, 'REQUESTED. TYPE YES (Y) TO CALCULATE, NO (N) 

OTHERWISE. ' 

READ (6, 99) REPLY 

IF (REPLY .EQ. 'Y' .OR. REPLY .EQ. 'y') S_FLAG = .TRUE. 

C SAMPLING TIME IS COMPUTED AS AVERAGE VALUE OF ALL DATA 

C THE COMPUTED VALUE CAN BE REPLACED BY OPERATOR - 

C IF OPERATOR SO DESIRES. 

C COMPUTE AVERAGE SAMPLE TIME FOR ALL DATA 

C DT= (LAST TIME - FIRST TIME)/ NUMBER OF POINTS MINUS 1 

C OR READ IN DT. 

2 DT = (TF - TO) /FLOAT (LF-1) 

PRINT * , ' DT IS : ' , DT , ' IS THIS OKAY TO USE - YES (Y) OR NO 

(N) ' 

READ (6, 99) REPLY 

IF (REPLY .EQ. 'N' .OR. REPLY .EQ. 'n') THEN 
PRINT *, 'ENTER DT IN SECONDS' 

READ* , DTNEW 
ELSE 

GO TO 4 
ENDIF 

C CHECK TO SEE IF THE NEW DT IS LEGAL 

IR = INT (0.1 + DTNEW/ DT) 

IF (IR.LT.2) THEN 

PRINT *, 'CANNOT REDUCE THE DT - MUST USE COMPUTED 

VALUE' 

GO TO 4 
ELSE 

1 = 1 
J = 1 

3 X ( I) = X(J) 

Y (I) = Y ( J) 

TLG(I) = TLG(J) 

IF ( J+IR .LT. LF) THEN 
1 = 1 + 1 
J = J + IR 
GO TO 3 
ELSE 

LF = I 
GO TO 2 
ENDIF 


3 


ENDIF 

4 CONTINUE 

TLNGTH = 0.5 *( TLG(l) + TLG(LF) ) 

CALL FBAND (FLOW, FHIGH, TLNGTH, PF_SRCH_BAND ) 

LEAST = INT ( 6./ (FLOW+FHIGH) ) 

C LEAST IS ESTIMATE OF HOW MANY POINTS TO USE FOR A MINIMUM 

C OF 3 CYCLES OF SKIPROPE. 

C 

c 

5 CONTINUE 

PRINT*, ' RECOMMEND USING AT LEAST ' , LEAST , ' POINTS FOR' 
PRINT*,' 3 DATA CYCLES. LAST POINT IN BUFFER IS ' , LF 
PRINT*, 'WHAT IS THE STARTING TIME INDEX ?' 

READ*, LB 

PRINT*, 'WHAT IS THE LAST TIME INDEX ? ' 

C 

READ* , LE 

IF (LE .GT. LF) THEN 

PRINT*, 'NOT ENOUGH POINTS IN BUFFER TO MAKE AN ACCURATE 

RUN. ' 

PRINT*, 'DO YOU WISH TO ENTER NEW START TIME AND LAST 

TIME' 

PRINT*, 'INDICES (Y OR N) ? IF NO, THE PROGRAM WILL ABORT 

AND' 

PRINT* , 'REQUEST REFILLING THE BUFFER AND RUNNING AGAIN.' 
READ (6, 99) REPLY 

IF (REPLY. EQ. 'Y' .OR. REPLY. EQ. 'y') GO TO 5 

PRINT*, 'PROGRAM WILL ABORT - REFILL BUFFER AND RUN AGAIN.' 
STOP 
END IF 
C 

IF ( 0 . EQ. MOD (LE-LB+1 , 2 ) ) LE=LE-1 

C THIS MAKES LE SUCH THAT NUMBER OF POINTS (LE-LB+1) IS ODD 

TLNGTH = 0.5 *( TLG(LB) + TLG(LE) ) 

7 PRINT *,' TETHER LENGTH IS :', TLNGTH,' KILOMETERS - OK TO 

USE' 

PRINT *,' REPLY WITH A YES (Y) OR A NO (N) ' 

READ (6,99) REPLY 

IF (REPLY .EQ. 'N' .OR. REPLY .EQ. 'n') THEN 
PRINT * , ' ENTER TETHER LENGTH IN KILOMETERS ' 

READ* , TLNGTH 
GO TO 7 
ENDIF 
C 

C CALL FBAND TO GET FREQUENCY SEARCH BANDS 

C 

CALL FBAND (FLOW, FHIGH, TLNGTH, P F_SRCH_B AND ) 


C 
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9 PRINT*, 'FLOW & FHIGH = ',FLOW, FHIGH 

PRINT*,' ARE THESE BOUNDS OKAY TO USE - Y OR N' 
READ (6,99) REPLY 

IF (REPLY .EQ. 'N' .OR. REPLY .EQ. 'n') THEN 
PRINT*, 'ENTER THE TWO VALUES IN HZ' 

READ*, FLOW, FHIGH 
GO TO 9 
END IF 

C 

99 FORMAT (Al) 


PRINT* ,' START INDEX - STOP INDEX - TOTAL POINTS PROCESSED ' 
PRINT*, LB, LE, LEB 

PRINT *,'DT = ' , DT , ' * * * TETHER LENGTH IN KILOMETERS = 

' , TLNGTH 

DF=1.0/ (NFT*DT) 

TSHIFT = DT* (LB-1) 

TMIDPT = TO + DT * ( (LB+LE-2)/2 ) 


C TSHIFT IS DELTA TIME FROM START TIME OF BUFFER (I.E. TO) 

C TO FIRST TIME POINT USED IN THIS RUN (DATA WINDOW) . 

Q 

C TMIDPT IS TIME POINT OF MIDDLE OF DATA WINDOW. 

C 


CALL WORK ( 11 , X, AMPX, PHASEX , FREQX, IWXMAX, G_FLAGX, AVGX, 

ODF FFT ) 

- CALL WORK ( 12, Y, AMP Y, PHASE Y , FREQY , IWYMAX , G_FLAGY, AVGY , 
ODF_FFT ) 

PSIGN =0.0 

C PSIGN SET TO ZERO - DEFAULT VALUE IN CASE FILTER CAN'T 

PREDICT 

C SKIPROPE 

C CHECK ON GOODNESS FLAGS 

C 

IF (G FLAGX . AND . G_FLAGY ) THEN 

FAVG = 0.5 * (FREQX + FREQY) 

TAVG =1.0/ FAVG 

WK = 1000. *TLNGTH*TAVG / (360*PI) 

UMAX = WK * AMPY 

VMAX = WK * AMPX , „ 

PRINT*, 'X AMP = ' , AMPX, ' X PHASE = PHASEX* 180 . 0./ PI , 

* ' X FREQ = ' , FREQX 

PRINT*, 'Y AMP = ' ,AMPY, ' Y PHASE = ' , PHASEY*180 . 0/PI , 
jf ' Y FREQ = ' , FREQY 

PRINT *, 'PHASE DIFFERENCE (DEGREES) BETWEEN X & Y = ' 

f , 57.3 * (PHASEX - PHASEY) 
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PRINT* , 'MAX U = UMAX,' MAX V = ', VMAX 
ELSEIF ( . NOT . G_FLAGY . AND . G_FLAGX) THEN 
FAVG = FREQX 
TAVG =1.0/ FAVG 

WK = 1000 . *TLNGTH*TAVG/ (360 * PI) 

VMAX = WK * AMPX 
UMAX = 7777. 


C 


PRINT * , ' Y AXIS DATA IS SUSPECT - FREQ OUT OF BAND' 

PRINT * 'FREQUENCY RETURNED IS PREDICTED MIDPOINT' 

PRINT * , ' DATA SHOULD NOT BE USED WITHOUT CAUTION' 

PRINT * , ' MAX V = ' , VMAX 

PRINT *, 'NEITHER POLARITY NOR YAW MANEUVER CALCULATIONS ARE 


PERFORMED ' 

GO TO 88 

ELSEIF ( . NOT . G_FLAGX . AND . G_FLAGY) THEN 
FAVG = FREQY 
TAVG =1.0/ FAVG 

WK = 1000. *TLNGTH*TAVG/ (360*PI) 

UMAX = WK * AMPY 


PRINT *^'X ? AXIS DATA IS SUSPECT - FREQ OUT OF BAND' 

PRINT * 'FREQUENCY RETURNED IS PREDICTED MIDPOINT' 

PRINT *,'DATA SHOULD NOT BE USED WITHOUT CAUTION' 

PRINT * 'MAX U = ' , UMAX 

PRINT * , ' NEITHER POLARITY NOR YAW MANEUVER CALCULATIONS ARE 


PERFORMED ' 

GO TO 88 

ELSEIF ( ( . NOT. G FLAGX) . AND. ( . NOT . G_FLAGY ) ) THEN 
PRINT * 'BOTH AXES ARE BAD ***** 3 SUGGESTIONS' 
PRINT *,'l) SUGGEST LOOK AT TIME PLOTS OF GYRO 


SIGNALS . ' 

PRINT *, '2) 
PRINT * , ' 3 ) 
UMAX = 7777. 
VMAX = 7777. 
FAVG = 7777. 
TAVG = 7777. 
WK = 7777. 
PSIGN= 7777. 


SUGGEST MAKE FFT PLOTS AND LOOK AT DATA.' 
SUGGEST WIDENING SEARCH BAND.' 


GO TO 88 


C 

C 

c 


DATA IS BAD. WRITE OUTPUT AT LABEL '88' 
END IF 


C 

C 

c 

c 

c 


PRINT *, 'AVERAGE PERIOD IN SECONDS IS :',TAVG 

THE INTEGERS 'IWXMAX' AND 'IWYMAX' ARE TIME INDICES WHERE X 
AND Y VALUES ARE A MAXIMUM. THIS CORRESPONDS TO WHERE 
COSINE (PHI) = 1 OR PHI = 2*PI. 

IF IWYMAX GT IWXMAX , MEANS POLARITY IS POSITIVE ABOUT Z 
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noon on oooo ojooo 


c 

c 

c 

c 

c 

c 


PROVIDED THE TIME DIFFERENCE BETWEEN IWYMAX AND IWXMAX IS 
EQUIVALENT TO 90 DEGREES. 

IWYMAX COULD BE GREATER THAN IWXMAX FOR NEGATIVE ROTATION 
BUT THIS WOULD REQUIRE A 270 DEGREE TRAVEL TIME. 

THUS THE TEST FOR POLARITY IS 180 DEGREES TRAVEL TIME. 

IF (IWYMAX .GT. IWXMAX) THEN 

TEST = DT * FLOAT (IWYMAX-IWXMAX) 

TEST IS TIME IN SECONDS TO GO FROM X-AXIS TO Y-AXIS. 
POLARITY DICTATED BY THIS TIME BEING GT OR LT 1/2 OF 


IF (TEST .GT. 0.5*TAVG) THEN 
PSIGN = -1.0 
ELSE 

PSIGN = +1.0 
ENDIF 
END IF 

NOW DO CASE FOR X PEAK OCCURS AFTER Y PEAK 
THIS IS SAME LOGIC AS ABOVE IN PRINCIPLE 

IF (IWXMAX .GT. IWYMAX) THEN 

TEST = DT*FLOAT( IWXMAX - IWYMAX) 

IF (TEST .GT. 0.5*TAVG) THEN 
PSIGN = +1.0 
ELSE 

PSIGN = -1.0 
ENDIF 

ENDIF 

ALL DONE - SIGN COMPUTATIONS ARE COMPLETED 
PRINT POLARITY OF SKIPROPE = PSIGN 

PRINT* 

WRITE TIME DATA TO FILE FOR RECORD ONLY IF REQUESTED. 
REQUEST IS IF ODF_TIME FLAG IS TRUE. 

IF (ODF_TIME) THEN 
DO I = 1, LEB 

TWX = COS ( 2 . 0*PI*FAVG*DT* (1-1) +PHASEX) 

TWY = COS ( 2 . 0*PI*FAVG*DT* (1-1) +PHASEY) 

WX = TWX * AMPX 
WY = TWY * AMPY 
U = -PSIGN * WK * AMPY * TWX 
V = -PSIGN * WK * AMPX * TWY 
T = TO + (LB + I - 2 ) *DT 
WRITE ( 13, *)T,WX,WY,U,V 

END DO 
ENDIF 

IF STATION 2 FLAG SET TO ' TRUE ' , THEN DO YAW MANEUVER 
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C CALCULATIONS, OTHERWISE SKIP TO LABEL 88 . 

C 

IF (S_FLAG) THEN 
C 

c CALCULATE NUMBER OF ROTATIONS ORBITER SHOULD EXECUTE. 

C 

RNROT = ( AMIN1 (UMAX, VMAX) )/ (2.0 * R_ARM) 

C 

c SPECIFY POLARITY 

IF (PSIGN .EQ. 1.) THEN 
POLARITY = 'POSITIVE' 

ELSE 

POLARITY = 'NEGATIVE' 

ENDIF 

C 

C OUTPUT REV NUMBER, POLARITY, # OF ROTATIONS, AND START 

TIMES 

C THIS DATA ALSO GOES TO FILE. DATA SET TIME TAG IS GIVEN 

BY 'TMIDPT' . 

C 

C BURN TIMES CALCULATED HERE ASSUME THAT THE ORBITER NOSE IS 

C ALIGNED WITH THE X-LVLH AXIS 

C IF THIS IS NOT THE CASE, THE BURN TIMES MUST BE ADJUSTED TO 

C ACCOUNT FOR WHERE THE NOSE IS WRT THE X_LVLH AXIS 

C THIS TIME ADJUSTMENT IS : (YAW ANGLE/360*FAVG) IN SECONDS. 

C YAW ANGLE IS DEGREES NOSE IS AWAY FROM X-LVLH. 

C FAVG IS VALUE FROM FREQUENCY MEASUREMENTS (IN HZ) . 

C 

C 

C CALCULATE 5 BURN TIMES AND OUTPUT TO SCREEN AND FILE 17 

C FIRST, ESTABLISH TIME WHEN X AND Y ARE MAXIMUMS - TXMAX 

& TYMAX 
C 

TXMAX = DT* (IWXMAX-1) 

TYMAX = TXMAX + 0.25*TAVG + TSHIFT + TO 
WRITE (6,76) TMIDPT, TF 
WRITE ( 6 ,*) 

WRITE (17,76) TMIDPT, TF 
WRITE (17,77) 

76 FORMAT (' MIDPOINT TIME FOR THIS DATA WINDOW IS 

• t # 5 f j 

$ ,' TIME TAG ON LAST POINT IN BUFFER IS :',E18.6) 

PRINT* 

77 FORMAT (' REV LABEL POLARITY RATE(D/S) # OF 
ROTATIONS ' 

$ ,' START TIME' ) 

C 

C NOW ADJUST BURN START TIMES TO ACCOUNT FOR ORBITER 

ORIENTATION. 

C THIS REQUIRES OPERATOR INPUT FOR ORBITER YAW ANGLE. 

C 

61 PRINT* ,' ENTER ORBITER YAW AXIS ANGLE IN DEGREES' 
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READ *, OYAWANG 

PRINT*, 'ORBITER NOSE IS OYAWANG, ' DEGREES WRT X-LVLH 

PRINT* , ' IS THIS CORRECT (Y OR N)?' 

READ (6, 99) REPLY 

IF (REPLY .EQ. 'N' .OR. REPLY .EQ. 'n') GO TO 61 


C 

C 

COMES 

C 

(TF) . 

C 

C 

78 


BTDEL = PSIGN * OYAWANG/ ( 3 60 . 0*FAVG) 

TYMAX = TYMAX + BTDEL 
WRITE (6,77) 

ADVANCE TIME 'TYMAX' BY INCREMENTS OF TAVG UNTIL 'TIME' 

UP THAT IS GREATER THAN TIME OF LAST DATA POINT IN BUFFER 

ANY COMPUTED TIMES LESS THAN 'TF' ARE IN THE PAST. 

IF (TYMAX .LE. TF) THEN 
TYMAX = TYMAX + TAVG 
GO TO 78 
END IF 


C 

C 

C 


17 


NOW COMPUTE THE 5 BURN TIMES 
DO K=1 , 5 

TT = TYMAX + (K-l ) *TAVG 

WRITE (6, 17) K-l, POLARITY, 360.*FAVG, RNROT, TT 
WRITE (17,17) K-l, POLARITY, 360.*FAVG, RNROT, TT 

END DO 

FORMAT ( 15 , 7X, A8 ,F12.4,7X,F5.1,7X, F12 . 2 ) 


END IF 


88 CONTINUE 


C 

C 

C 

C 

C 


CALL ODTF (TO, TMIDPT, TF, DT, LE, LB, LEB, TLNGTH 
$ , AMPX, FREQX, PHASEX, AMPY, FREQY , PHASEY 

$ , FAVG, TAVG, WK, PSIGN, UMAX, VMAX, FLOW, FHIGH 

$ , AVGX, AVGY) 


CLOSE ALL FILES 


CLOSE 

CLOSE 

CLOSE 

CLOSE 

CLOSE 

END 


(11, STATUS= ' KEEP ' ) 
(12, STATUS= ' KEEP ' ) 
(13, STATUS= ' KEEP ' ) 
(17, STATUS= ' KEEP ' ) 
(18, STATUS='KEEP' ) 
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C SUBROUTINE WORK CALCULATES THE AMPLITUDE, PHASE , AND 

FREQUENCY 

C OF THE DATA. THE FOURIER TRANSFORM SUBROUTINE FOUR1 IS 

C CALLED BY SUBROUTINE WORK. WORK RETURNS TO THE MAIN 

PROGRAM 

C THE VALUES OF THE AMPLITUDE , PHASE, AND FREQUENCY AS WELL AS 

C THE TIME INDEX WHERE THE MAXIMUM VALUE OCCURS. 

C THIS IS BASED ON MODEL OF COS (WT+PHASE) . 

C 

SUBROUTINE WORK (IOA, ANG, AMP, PHASE, FREQ, ITMAX 
$, G_FLAG , BIAS, FFT_FLAG ) 

INTEGER NDIM, NCDIM 

PARAMETER (NDIM=3000, NCDIM=8200, NFT=8192, NPNT=7) 
PARAMETER (PI=3 . 14 15926 , DFR=57 . 2957795 , RFD=0 . 017453292 ) 
DIMENSION AUX (NDIM) , ANG (1) 

REAL *4 XFREQ(7) , PHIMAG ( 7 ) , PHREAL ( 7) 

COMPLEX AWO (NCDIM) 

LOGICAL G FLAG , FFT_FLAG 


COMMON LB, LE, DT, DF, FLOW, FHIGH 
NTB1=LE-LB+1 

C NTB1 IS FORCED TO BE ODD IN MAIN PROGRAM. 

C HANN WINDOW ROUTINE USES ODD NUMBER OF POINTS TO TAPER. 

C 

C 

C LOAD INPUT DATA FROM ANG (I) INTO ARRAY AUX(J). 

LB1=1-LB 
DO I=LB, LE 
IL=I+LB1 
AUX ( IL) =ANG ( I ) 

END DO 

C 

C APPLY WINDOW FUNCTION TO TIME SEQUENCE 

CALL HANN (NTB1 , AUX , BIAS) 

C 

C MAKE COMPLEX NUMBER AWO (I) FROM REAL NUMBER AUX (I) BY USING 

C A ZERO IMAGINARY VALUE (AUX (I) IS THE REAL VALUE). 

C 

DO 1=1 , NTB1 

AWO ( I ) =CMPLX (AUX ( I ) ,0. ) 

END DO 

C 

C NOW PAD THE DATA STREAM WITH ZEROS OUT TO AWO (8192) . 

C 
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DO I=NTB1+1 , NFT 

AWO(I) = CMPLX ( 0 . ,0. ) 
END DO 


C 

C SUBROUTINE FOUR1 DOES THE FOURIER TRANSFORM USING A FTT 

METHOD . 

C 

CALL FOUR1 (AWO,NFT, 1) 

C NOW FIND THE MAXIMUM MODULUS VALUE OF THE FOURIER TRANSFORM 

DATA 

C OVER A SPECIFIED FREQUENCY INTERVAL. THIS INTERVAL IS 

CALCULATED 

C FROM INPUT DATA AND IS DESIGNED SUCH THAT THE SKIPROPE 

C FREQUENCY FALLS WITHIN THIS INTERVAL. THE INTERVAL IS 

C SUFFICIENTLY NARROW THAT NO OTHER MODE SHOULD FALL WITHIN 

THE 

C INTERVAL. 


C FLOW IS LOWER BOUNDARY OF SEARCH BAND 

C FHIGH IS UPPER BOUNDARY OF SEARCH BAND 

C 

XMAX=0 . 0 

IFRST = 1 + I NT ( FLOW/DF) 

I LAST = 1 + INT ( FHIGH/DF) 

DO I = IFRST, I LAST 
FR = (1-1) *DF 

IF (CABS (AWO (I) ) .GT.XMAX) THEN 
XMAX=CABS (AWO (I) ) 

KF=I 

FREQ = FR 
END IF 
END DO 


C CREATE THE 3 DATA SETS TO BE FITTED BY LEAST SQUARES 

POLYNOMIAL . 

C POLYNOMIAL IS 2ND DEGREE AND 7 POINTS WILL BE USED IN CURVE 

FIT. 

C 3 SETS ARE: 

C MAGNITUDE OF TRANSFORM ( SQRT ( REAL* * 2 + IMAG**2) ) 

C REAL PART 

C IMAGINARY PART 

C 

C CENTER OF DATA SET IS THE FREQUENCY POINT WHERE MAX WAS 

FOUND. 

C 

DO I = 1 , NPNT 

J = KF- ( (NPNT+1) /2 . 0) +1 
XFREQ(I) = CABS ( AWO ( J) ) 
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PHIMAG(I) = AIMAG (AWO ( J) ) 
PHREAL(I) = REAL (AWO (J) ) 
END DO 


C 

C 

c 

c 


DO CURVE FIT ON THE MODULUS OF THE FOURIER TRANSFORM 
CALL TO LSCF WITH OPTION 1 DOES 2 THINGS. 

CURVE FITS AND COMPUTES TRUE MAX FREQUENCY POINT. 


CALL LSCF ( FQ_P0 , XMAX, XFREQ, 1, G_FLAG) 

C G_FLAG = TRUE MEANS MAX FREQUENCY FOUND IN THE SPECIFIED 

INTERVAL^ is TH£N CALLED twice with OPTION 2 TO EVALUATE THE 

polynomial^ critical frequency value found IN THE FIRST LSCF 

C ALL ' IF FALSE, THEN THE MAXIMUM PEAK OCCURS OUTSIDE THE 7 POINT 
^ ANGE ‘ (THIS IS POSSIBLE IF THE INITIAL ESTIMATE OF THE SKIPROPE 
FREQUENCY ^ TRTRER LENGTH is MUCH DIFFERENT FROM THE 

TRUE VALUE.) WHEN G_FLAG IS FALSE THE INFORMATION RETURNED 
VALUES BASED ON THE MIDPOINT OF THE SPECIFIED SEARCH BAND. 


C 

ARE 

C 


IF (G FLAG) THEN 

CALL LSCF (FQ_PO, PHASEI, PHIMAG, 2, G_FLAG) 
CALL LSCF ( FQ_PO , PHASER, PHREAL, 2, G_FLAG) 

EL SE . . 

KF = 1+ INT( ( (FLOW+FHIGH) /2.0 ) /DF ) 

XMAX = CABS (AWO (KF) ) 

PHASEI - AIMAG ( AWO(KF)) 

PHASER = REAL( AWO(KF)) 

FREQ = 0 . 

FQ_PO = KF - 1.0 
END IF 

FREQ = FREQ + DF * FQ_P0 


C SCALING OF TRANSFORMED DATA IS PERFORMED TO GIVE OUTPUTS IN 

C DEGS/SEC AND REPRESENT ACTUAL RATE DATA. 

SCALE = 4 . 0/ FLOAT (NTB1-1) 

AMP = SCALE * XMAX 

PHASE = -ATAN2 (PHASEI, PHASER) 

C THE FORWARD FOURIER TRANSFORM IS USUALLY DEFINED WITH 

C XP(_;L MANY F FFt‘rOUTINES, INCLUDING FOUR1 , USE EXP (+i*2*PI*F*T) . 

THESE ^ rferent CONVENTIONS FOR THE FORWARD FOURIER TRANSFORM 
RESULT IN TWO 
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c DIFFERENT FORMS FOR THE SHIFT THEOREM. IN THE FIRST CASE, 

THE SHIFT 

c THEOREM STATES THAT IF G (T) TRANSFORMS AS G(F), THEN 

G (T+Tl) TRANSFORMS „ _ 

q EXP (i*2*PI*F*Tl) *G(F) . IN THE SECOND CASE, IF G(T) 

TRANSFORMS AS . 

c g (F) , THEN G (T+Tl) TRANSFORMS AS EXP (-1*2*PI*F*T1) *G (F) . 

SINCE O UR 

c MODEL IS COS ( 2 *PI*F*T + P) = COS (2*PI*F* (T + P/ (2*PI*F) ) , 

AND WE USE 

C THE FIRST CONVENTION FOR THE FOURIER TRANSFORM, WE EXPECT 

OUR PHASE 

c TO BE 2 *PI*F*P/ ( 2 *PI*F) = P. HOWEVER, SINCE THE PROGRAM 

USES THE 

c SECOND CONVENTION FOR THE FOURIER TRANSFORM, THE PHASE IS 

— P SO TO 

c ' CORRECT FOR THIS DIFFERENCE WE MUST INCLUDE ANOTHER - SIGN: 
-(-P) = P. 


CALCULATE THE TIME INDEX WHERE THE MAXIMUM RATE OCCURS. 
THIS IS BASED ON COS (PHI) = 1 IMPLIES PHI = 2*PI 

ITMAX=INT ( (1.0 / (FREQ*DT) ) * ( 1 . 0-PHASE/ (2 . 0*PI) ) +0.5) +1' 

IF THIS INDEX CORRESPONDS TO A TIME GREATER THAN 1 PERIOD, 
THEN SUBSTRACT THE EQUIVALENT OF 1 PERIOD FROM ITMAX. 

IF ( ITMAX*DT .GT. (l./FREQ) ) THEN 

ITMAX ■ ITMAX - INT ( 1.0/ (FREQ*DT) ) 

ENDIF 


C OUTPUT THE MODULUS OF THE TRANSFORM FROM 0.0 HZ THROUGH THE 

PENDULOUS 

C FREQUENCY (ASSUMING MAX VALUE IS 0.04 HZ), USING A SPACING 

OF DF ♦ 

c OUTPUT 4 NUMBERS PER LINE: FREQ (HZ) , MODULUS (DEG/ SEC) , REAL 

PART 

c ' AND IMAGINARY PART. (NOTE: LAST TWO ARE NOT IN DEG/ SEC) 
DF=1. / (8192*8*0.128) = 1./8388.6 = 0.0001192 
KQ1 = . 04 /DF = 335.54 CALL THIS 336 


C 

C 

C 

C 

c 

c 


WRITE FFT DATA TO OUTPUT FILE IF FFT_FLAG IS TRUE 
OTHERWISE DO NOT WRITE TO OUTPUT. 


IF (FFT_FLAG) THEN 
I LAST = 1006 

IF (FREQ .GT. .0035 ) I LAST = 336 
DO I = 1, I LAST 
FR = (I — 1 ) *DF 
XMAX = SCALE*CABS (AWO (I) ) 

WRITE (IOA,*) FR, XMAX, REAL(AWO(I) ) , AIMAG(AWO(I) ) 
END DO 
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ENDIF 

RETURN 

END 


SUBROUTINE HANN (LA, All, BIAS) 


TAPER IS RAISED COSINE CURVE. 

MEAN IS COMPUTED AND REMOVED FROM INPUT SIGNAL 
PARAMETER (PI=3 . 1415926) 

REAL All (LA) , HW(3000) 

ITM = (LA-1) / 2 
RM = FLOAT (ITM) 

DO IT= -ITM, ITM 

I = 1 + IT + ITM 

HW(I)=0.5*(1.0 + COS (PI*FLOAT (IT) /RM ) ) 

A11(I)=A11(I) *HW(I) 


END DO 

COMPUTE MEAN 
TRUE MEAN IS 
REDUCES VALUE 


OF TAPERED SIGNAL 

TWICE COMPUTED VALUE BECAUSE HANN WINDOW 
BY FACTOR OF 2. (I.E. MEAN OF WINDOW IS 


0.5) 


CALL MEAN (LA, All, BIAS) 

BIAS = BIAS * 2.0 * LA/ (LA-1) 

DO I = 1,LA 

All ( I ) = All (I) ~ BIAS * HW (I) 

END DO 
RETURN 
END 


SUBROUTINE MEAN (LA, A22 , SA) 

THIS ROUTINE COMPUTES THE DC TERM OF THE DATA STREAM. 
MEAN IS NOT REMOVED, BUT ONLY COMPUTED. 

REAL A22 (LA) 

SA = 0. 

DO 1=1, LA 

SA=SA+A22 (I) 

END DO 

SA=SA/ FLOAT (LA) 

RETURN 

END 


SUBROUTINE FOUR1 (DATA, NN, ISIGN) 

THIS ROUTINE DOES THE FOURIER TRANSFORM USING A FFT METHOD. 

REAL* 8 WR , WI , WPR , WPI , WTEMP , THETA 
DIMENSION DATA ( * ) 
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N=2*NN 

J=1 

DO 11 1=1, N, 2 
IF(J.GT.I)THEN 
TEMPR=DATA(J) 

TEMP I=DATA( J+l) 

DATA ( J ) =DATA ( I ) 

DATA(J+1)=DATA(I+1) 

DATA ( I ) =TEMPR 
DATA(I+1)=TEMPI 
ENDIF 
M=N / 2 

1 IF ( (M.GE.2) .AND. (J.GT.M) ) THEN 

J=J-M 
M=M/2 
GO TO 1 
ENDIF 
J=J+M 

11 CONTINUE 

MMAX=2 

2 IF (N.GT.MMAX) THEN 

ISTEP=2*MMAX 

THETA=6. 283 185307 17959D0/ (ISIGN*MMAX) 

WPR=-2 .D0*DSIN(0.5D0*THETA)**2 
WPI=DSIN (THETA) 

WR=1 . DO 
WI=0 . DO 

DO 13 M=1 ,MMAX, 2 
DO 12 I=M, N , I STEP 

TEMPR=SNGL(WR) *DATA( J) -SNGL(WI) *DATA(J+1) 

TEMPI=SNGL(WR) *DATA( J+l) +SNGL(WI) *DATA( J) 

DATA ( J ) =D AT A ( I ) -T EMPR 

DATA ( J+ 1 ) =DATA ( 1+1 ) -TEMPI 

DATA ( I ) =DATA ( I ) +TEMPR 

DATA ( 1+1) =DATA ( 1+1) +TEMPI 

12 CONTINUE 
WTEMP=WR 

WR=WR*WPR-WI*WPI+WR 
WI=WI *WPR+WTEMP* WPI+WI 

13 CONTINUE 
MMAX=ISTEP 

GO TO 2 
ENDIF 
RETURN 
END 


C THIS SUBROUTINE DOES LEAST SQUARES CURVE FIT TO 7 POINTS 

C FOR A 2ND DEGREE POLYNOMIAL. THE DATA IS ASSUMED TO BE 

SAMPLED^ integral INTERVALS. ANY SCALING MUST BE DONE OUTSIDE 
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c 

c 

c 

c 

c 


THIS SUBROUTINE. THE 7 POINTS ARE : 

P = —3, “2, “1/ 0, 1, 2, 3 
THE POLYNOMIAL IS F(P) = A + B*P + C*P*P. 

THE MAX OCCURS AT P = PO = -B / (2*C) . 

FREQUENCY CORRESPONDING TO PO IS P0*DF (DF OF DATA 


STREAM) 

C 

POINTS. 


THIS DELTA IS REFERENCED TO MIDPOINT FREQUENCY OF 7 
THE MAX VALUE IS F(P0) = A - (B*B)/(4*C). 


SUBROUTINE LSCF (PO, FMAX, U_IN, IOPT, G_FLAG) 

ON ENTRY: 

U IN IS INPUT ORDINATE VALUES. (7 ) 

IOPT IS OPTION FOR 1 OF 2 THINGS 

1 : FIND PO WHERE MAX OCCURS PLUS COMPUTE MAX VALUE. 

2 : COMPUTE VALUE OF POLYNOMIAL AT SPECIFIED FREQUNCY 


C 

c 
c 
c 
c 
c 

PO. 

c 

POLYNOMIAL. 
C 


IF I0PT=2 , THEN PO IS FREQUENCY POINT TO EVALUATE 


C 

C 

C 

C 

C 

C 

C 

C 


C 

c 

c 

c 

c 

c 


ON EXIT „ 

PO IS VALUE OF P WHERE MAX PEAK OCCURS; 

THIS IS WRT CENTER POINT OF DATA. 

FMAX IS VALUE OF FUNCTION AT P=P0. 

G FLAG IS FLAG FOR DATA VALIDITY 
“ SET TO TRUE IF EVERYTHING IS OKAY 

SET TO FALSE IF PEAK IS OUTSIDE SEARCH ZONE 
******************************************** 

REAL* 4 U_IN ( * ) 

LOGICAL G FLAG 

******************************************** 

FIRST STEP IS TO DO LEAST SQUARES. 

ALL COEFFICIENTS HAVE BEEN PRE-COMPUTED . 

'A' IS -8, 12, 24, 28, 24, 12, -8 DIVIDED BY 84 

'B' IS -9, -6, -3, 0, 3, 6, 9 DIVIDED BY 84 

/ C / IS 5, 0, -3, -4, -3, 0, 5 DIVIDED BY 84 

US17 = U_IN ( 1 ) + U_IN (7 ) 

US35 = U_IN ( 3 ) + U_IN(5) 

A = COF1 

COF1 = -8 . *US17 + 12 . * (U_IN(2) +U_IN(6) ) + 

^ J -if- . ^ r» T7 TXT t A 


# 


B = COF2 
COF2 = 9 


# 

* 


C= COF3 
COF3 = 5 


24 . *US35 + 28 . *U_IN (4) 

*(-U IN (1) +U_IN ( 7 ) ) + 

6 . * ( -U_IN ( 2 ) +U_IN ( 6 ) ) + 

3 . * ( -U_IN ( 3 ) +U_IN ( 5 ) ) 


*US17 -3 . *US35 - 4 . *U_IN (4 ) 

IF (ABS (COF3) .LT. 1. 0E-08) THEN 

PRINT* / ********* WARNING ********** 

PRINT*' 'CONSTANT VALUE EQUALS ZERO - CANNOT COMPUTE A 
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MAXIMUM' 

PRINT* FREQUENCY VALUE.' 

G_FLAG = .FALSE. 

RETURN 
ENDIF 

DID NOT DIVIDE BY 84 YET. DO SO FOR MAX PART BUT NOT PO . 


COMPUTE PO, VALUE OF F(P) WHERE A+B*P+C*P*P - 0 
COMPUTE FUNCTION AT PO; A+B*PO+C*PO*PO = A-B**2/4C 
IF ( IOPT .EQ. 1) THEN 
P0=-0 . 5*COF2/COF3 

FMAX = (COF1 + 0.5 * PO * COF2) / 84. 

G FLAG = .TRUE. 


IF ( ABS (PO) .GT. 3.0) THEN 

PRINT*,'* * * WARNING * * *' 

PRINT*, 'HAVE NO MAX VALUE IN SPECIFIED 
G_FLAG = .FALSE. 

ENDIF 

ELSE „ . 

FMAX = (COF1 + PO* ( COF2 + PO*COF3 ) )/84. 

END IF 
RETURN 
END 


INTERVAL' 


SUBROUTINE FBAND (FL, FH, TLKM, PF) 

C SUBROUTINE RETURNS FL AND FH; THESE ARE LOW AND HIGH 

FREQUENCY 

n vattifs TO SEARCH FOR PEAK AMPLITUDE. 

C THIS SUBROUTINE COMPUTES SKIPROPE FREQUENCY .(FC) FROM THE 

r SES OF THE ORBITER AND SATELLITE AND THE TETHER LENGTH. DATA 

NEEDED^^ THESES CALCULATIONS ARE IN COMMON BLOCK 'FREQ' AND ARE 
READ 

c FROM FILE 'PARAM.DAT' IN MAIN PROGRAM. 

C TKLM IS THE TETHER LENGTH IN KILOMETERS. ™ 

C PF IS PERCENT OF SKIPROPE FREQUENCY TO USE AS A DELTA F , 

C* E *' BAND IS FROM FC - DELTA (FL) TO FC + DELTA (FH) . 

p 

REAL* 4 DENSITY, TOTALL, MSAT , MORB, ALTKM, FL, FH, FC, 

TLKM PF 

COMMON/ FREQ/ DENSITY, TOTALL, MSAT, MORB, ALTKM 
OMSQ = ORBRATESQ (ALTKM) 

CK = 3.0 * OMSQ / DENSITY 

MO - MORB + TOTALL * DENSITY 

Q = 0.5 * DENSITY * TLKM 

MSTAR = ( (MO-Q) * (MSAT+Q) ) / (MO+MSAT) 
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FC = 0.5 * SQRT (CK*MSTAR/TLKM) 


DF = FC*PF / 100 . 

FL = FC - DF 
FH = FC + DF 

RETURN 

END 

REAL FUNCTION ORBRATESQ (ALTKM) 
PARAMETER (GM = 9.81098, RE = 6378.17) 
R = GM/ (1000. 0*RE) 

ORBRATESQ = R/(1.0 + ALTKM/ RE) **3 

RETURN 

END 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

2 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE READINDATA (TO, TF , X, Y, TLG, MODE, LF, MFLAG) 

SUBROUTINE READS IN DATA FROM FILE. THIS EQUIVALENT 
TO DATA THAT WILL COME FROM THE PREPROCESSOR BLOCK. 


REAL* 4 X ( 1) , Y ( 1) , TLG ( 1) 
INTEGER MODE 
LOGICAL MFLAG 


TO 

TF 

X(I) 

Y ( I ) 
TLG (I) 
MODE 


LF 

MFLAG 


' TF ' 


?IME TAG FOR 1ST POINT IN BUFFER 
?IME TAG FOR LAST POINT IN BUFFER 
C AXIS GYRO DATA (DEG/SEC) 
r AXIS GYRO DATA (DEG/SEC) 

■“ETHER LENGTH IN KILOMETERS 
U4CSMODE VALUE FOR LAST TIME POINT 
AMCSMODE =0, 1, 2, OR 3 

0 - NO VALID DATA 

1 - PASSIVE CASE 

2 - YAW HOLD 

3 - SPIN CASE 

DUMBER OF TIME POINTS STORED IN X & Y ARRAYS 
LOGICAL TO DENOTE IF AMCSMODE CHANGED FROM 1 OR 

TO A 0 OR 3 BETWEEN TIMES OF TO TO TF. 


AMCSMODE = 
AMCSMODE = 
AMCSMODE = 
AMCSMODE = 


JNOMSC . FOR READS FROM FORTRAN FILE 10, 5 NUMBERS 
TIME X-GYRO RATE, Y_GYRO RATE, TETHER LNGTH , 
TIME DATA SHOULD BE IN SECONDS. 

GYRO RATE DATA SHOULD BE IN DEGS/SEC. 
TETHER LENGTH IN KILOMETERS. 

MODE IS INTEGER BETWEEN 0 AND 3. 


PER LINE. 
AMCSMODE 


MFLAG = . TRUE . 

MFLAG IS LOGICAL VARIABLE TO DENOTE THAT MODE STATUS IS 
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c CONSISTENT FOR THE FULL DATA STREAM. IF MODE CHANGES TO 

C 0 (NO VALID DATA) OR 3 (SPIN) , THEN THE READ OPERATION IS 

C STOPPED AND THE BUFFER WILL HAVE A REDUCED NUMBER OF 

POINTS^e OPERATOR IS NOTIFIED BY A MESSAGE PRINTED TO THE 
SCREEN . 

C OPEN ( 10, FILE='IDFTXY. DAT ' ,STATUS=' UNKNOWN') 

C READ FIRST LINE TO GET FIRST TIME 

READ (10, *, END=2 ) TO, X(l), Y(l), TLG(l) , JMODE 
C LOOP TO READ DATA 

1 I READ ( 10 , * , END=2 ) TF, X(I), Y(I), TLG(I) , MODE 

1 = 1 + 1 

IF (JMODE .EQ. MODE) THEN 
GO TO 1 
ELSE 

IF ( (MODE* JMODE) .EQ. 2) THEN 
JMODE = MODE 
GO TO 1 
ELSE 

MFLAG = .FALSE. 

GO TO 2 
ENDIF 
ENDIF 

2 LF=I-1 

CLOSE (10, STATUS= ' KEEP ' ) 


RETURN 

END 


SUBROUTINE ODTF (TO, TM, TF , DT, LE, LB, LEB, TL 
$ , AX, FX, PX, AY, FY, PY 

$ , FA, TA, WK, PSIGN, UMAX, VMAX, FL, FH 

$ , A VGX, AVGY) 

WRITE DATA IN ARGUMENT TO FILE 

DEFINITION OF DATA ELEMENTS 

1ST TIME POINT IN BUFFER (MET) 

TIME AT MIDPOINT OF DATA WINDOW (MET) 

LAST TIME POINT IN BUFFER (MET) 

SAMPLE TIME - SECONDS 
INDEX OF LAST POINT IN DATA WINDOW 
INDEX OF FIRST POINT IN DATA WINDOW 
NUMBER OF POINTS USED IN DATA WINDOW 
TETHER LENGTH FOR THIS CASE - KILOMETERS 
PEAK MAGNITUDE OF RATE IN X AXIS - DEG/ SEC 


C TO 

C TM 

C TF 

C DT 

C LE 

C LB 

C LEB 

C TL 

C AX 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


FX 

PX 

AY 

FY 

PY 

FA 

TA 

WK 

PSIGN 

UMAX 

VMAX 

FL 

FH 

AVGX 

AVGY 


MEASURED FREQUENCY OF SKIPROPE IN X AXIS - HZ 
PHASE ANGLE IN X AXIS 

PEAK MAGNITUDE OF RATE IN Y AXIS - DEG/ SEC 
MEASURED FREQUENCY OF SKIPROPE IN Y AXIS - HZ 
PHASE ANGLE IN Y AXIS 
AVERAGE FREQUENY OF 'GOOD AXES' 

AVERAGE PERIOD OF 'GOOD AXES' 

CONVERSION CONSTANT FROM RATE TO U ; METER-SEC/DEG 
DIRECTION OF SKIPROPE - WRT Z LVLH AXIS 
MAXIMUM VALUE OF MIDNODE IN X DIRECTION 
MAXIMUM VALUE OF MIDNODE IN Y DIRECTION 
LOWER BOUNDARY OF FREQUENCY SEARCH BAND - HZ 
UPPER BOUNDARY OF FREQUENCY SEARCH BAND - HZ 
COMPUTED MEAN OF X AXIS TIME SEQUENCE 
COMPUTED MEAN OF Y AXIS TIME SEQUENCE 


WRITE (18,*) TM 

WRITE (18,*) LB, LE, LEB 

WRITE (18,*) TO, TF, DT, TL, WK, PSIGN 

WRITE (18,*) AX, FX, PX, AY, FY, PY 

WRITE (18,*) FA, TA, UMAX, VMAX, FL, FH 

WRITE (18,*) AVGX, AVGY 

RETURN 

END 
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APPENDIX 2 


APPENDIX 2. A 


Model Signal Generation Programs with Random Number Generator and 
Box-Muller Algorithm for Gaussian Distributed Variates 


APPENDIX 2 • B 


Results of ECR Verification Table 2. Filter Test Cases 

APPENDIX 2 . C 


Programs for ECR Testing . 

/ 1 } SIMREAD - Read Simulation Data Files . • 

(2) SIMTEST - Compare Observer Output to Simulation Inpu 

(3) BTBUSO - Compare Observer Output to Model Signal Input 


APPENDIX 2 . D 


Programs 

( 1 ) 

( 2 ) 

(3) 

(4) 


for Systematic Testing 

NO NOISE - Noise-free Test Cases, Station 2 
NOISE - Noisy Test Cases, Station 2 
LIBRATION - Noise-free Tests at. Station 1 
LIB NOISE - Noisy Tests at Station 1 


APPENDIX 2 . E 


Simulation Test Results 


APPENDIX 2.F 


Systematic Test Results ... , 

(1) NO NOISE - Noise-free Test Cases, Station 2 

(2) NOISE - Noisy Test Cases, Station 2 

(3) LIBRATION - Noise-free Tests at. Station 1 

(4) LIB_NOISE - Noisy Tests at Station 1 



APPENDIX 2. A 


Model Signal Generation Programs with Random Number Generator and 
Box-Muller Algorithm for Gaussian Distributed Variates 



The following random number generator is documented 
on page 196 of "Numerical Recipes" (Press, Flannery, Teukolsky, 
and Vetterling, 1989, Cambridge University Press). Setting th 
init ia 1 va lue idum to a negative integer initializes the routine. 


function rani (idum) 
dimension r(97) 

parameter (ml = 259200 , ial — 7141, 

parameter (m2 = 134456, ia2 = 8121, 

parameter (m3 = 243000, ia3 = 4561, 

data iff /0/ 

if (idum. lt.O.or. iff .eq.0) then 


icl = 54773, rml = 1.0/ml) 
ic2 = 28411, rm2 = 1.0/m2) 
ic3 = 51349) 


iff = 1 

ixl = mod (icl - idum, ml) 

ixl = mod(ial*ixl + icl, ml) 

ix2 = mod (ixl, m2) 

ixl = mod ( ial*ixl + icl, ml) 

ix3 = mod (ixl, m3) 

do j = 1,97 

ixl = mod(ial*ixl + icl,ml) 
ix2 = mod(ia2*ix2 + ic2,m2) 
r ( j ) = (float (ixl) + float (ix2) *rm2) *rml 


end do 
idum = 1 
end if 

ixl = mod(ial*ixl + icl, ml) 
ix2 = mod(ia2*ix2 + ic2,m2) 
ix3 = mod(ia3*ix3 + ic3,m3) 
j = 1 + (97*ix3) /m3 
if ( j .gt. 97 . or . j . It. 1) pause 

rani = r ( j ) . . , 

r ( j ) = (float ( ixl) + float (ix2) *rm2) *rml 


return 

end 


The following Box-Muller algorithm for generating 
normally (Gaussian) distributed variates is documented on page 
202 of "Numerical Recipes". The random number generator r 
listed above is used in this algorithm. 


function gasdev(idum) 
data iset/0/ 
if (iset.eq. 0) then 

vl = 2 . 0 *ranl (idum) - 1.0 
v2 = 2 . 0 *ranl(idum) - 1.0 
r = vl**2 + v2**2 
if (r.ge.l) go to 1 


1 


gset = vl*f ac 
gasdev = v2*fac 
iset = 1 

else 

gasdev = gset 
iset = 0 
end if 
return 
end 


The following two programs, CREATE. FOR and 
CRELIBR.FOR, use the random generator rani and the Box- Muller 
algorithm to generate Gaussian distributed noise to be added to 
the model test signals for the frequency domain skiprope 
observer, if noise is desired. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Program CREATE. FOR generates model data sets for the UNOMSC. FOR backup skip 
rope observer program. Five numbers are outputted per line of the data file 
IDFTXY.DAT: time, x (roll) axis gyro rate, y (pitch) axis gyro r ^ te ' 
length, and amcsmode. The user can opt to include noise and/or dropouts into 
the data If desired, the program will estimate the skiprope and pendulous 
frequencies as functions of the tether length. (Note: These estimated fre- 
quencies are exactly that - estimates - and definitely should not be con 
strued as accurate values!) 

The operator is prompted to input the following values: 
time sampling interval, beginning time tag 

number of data points, standard deviation, and random seed 

(standard deviation is 10 X the nominal value of 2.8e-04) 
(random seed should be a negative seed) 

skiprope amplitude, frequency, and phase (for x axis) 

skiprope amplitude, frequency, and phase (for y axis) 

1 or 2 to add noise or not, dropout percentage (0.0 to 1.0) 
orb rate, pendulous amplitude, frequency, and phase (for x axis) 

orb rate, pendulous amplitude, frequency, and phase (for y axis) 

tether length, amcsmode 

0 to use inputted values to generate model, 1 to estimate skiprope 
and pendulous frequencies from the inputted tether length 
If noise is added to the data, the program prints to the operator the 
values of the signal-to-noise ratios on each axis, XSNR and YSNR. 

The signal-to-noise ratio is calculated as the maximum value divided 
by the root mean square value for that axis. The noise is calculated 
using the Box-Muller method of generating Gaussian noise, and a 
portable random number generator called rani found in "Numerical 
Recipes" . 


C OPEN OUTPUT FILES 'IDFTXY.DAT' (FOR UNOMSC. FOR) AND 'F_REC.DAT' (FOR THE 
C ERROR CALCULATION PROGRAM VTBUSO.FOR) 

OPEN (9, FILE = 'F_REC.DAT' , STATUS = 'UNKNOWN') 

OPEN (10, FILE = 'IDFTXY.DAT' , STATUS = 'UNKNOWN') 

C READ INPUTS 
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PRINT* , ' ENTER THE TIME SAMPLING INTERVAL DT AND BEGINNING TIME 

PRINT^ENTER NUMBER OF TIME POINTS, STANDARD DEVIATION, RANDOM SEED' 
PRINT** ' (RANDOM SEED SHOULD BE A NEGATIVE NUMBER) ' 

PMNT^'SKIPROPE E aSpLITUDE ?, FREQUENCY ?, PHASE ? (FOR X)' 

PRINT*^' SKIPROP^AMPLITUDE ? , FREQUENCY?, PHASE? ( FOR Y ) 

PRINT^' NOISE* ? H YES « !, NO - 2; % DROPOUT DESIRED (0.0 TO 1.0)' 

PRINT*^ ENTER^ORB RATE, PENDULOUS AMP, FREQ, AND PHASE (FOR X)' 

PRINT*^ ENTER * ORB * RATE , PENDULOUS AMP, FREQ, AND PHASE (FOR Y) ' 

PRINT* , 'ENTER* THE TETHER LENGTH IN KILOMETERS AND AMCSMODE' 

READ* , TLNGTH , MODE 


C WRITE THE OUTPUT FILE 'F_REC.DAT' 

WRITE (9 , *) NOISE, SD, ISEED 
WRITE (9 , *) A1X , FIX, PHI1X 
WRITE ( 9 , * ) A1Y , F1Y , PHI1Y 

WRITE (9 , *) A2X, F2X, PHI2X 
C CONVERT PHASES FROM RADIANS TO DEGREES 
PI = 4 . 0*ATAN (1 . 0) 

CONVRT = PI/180.0 
PHI IX = PHI1X*C0NVRT 
PHI2X = PHI2X*CONVRT 
PHI1Y = PHI1Y*C0NVRT 
PHI2Y = PHI2Y*CONVRT 
TPIDT = 2 . 0*PI*DT 


C ESTIMATE THE SKIPROPE AND PENDULOUS FREQUENCIES FROM THE TETHER LENGTH 


0.2119024 
-0.3571371 
0.5309507 

SKIPC2**2 - 4 . 0*SKIPC3* (SKIPC1-TLNGTH) 

40. 0* ( (-1) *SKIPC2 + SQRT (DISCR1) ) / SKIPC3 
8 . 0952965E-02 
0.5571405 

0.2476189 

PENDC2**2 — 4 . 0*PENDC3 * (PENDC1 -TLNGTH) 

0.01*(1.0+((-l)*PENDC2 + SQRT (DISCR2) ) / ( 2 . 0*PENDC3 ) ) 
(TLNGTH. LT. 1. 2) THEN 

PENDFR = PENDFR + 0.008 
SKIPPR = 116. 049*TLNGTH 
END IF 

PRINT*, 'EST.^ SKIPROPE FREQ = ',SKIPFR,' EST. PENDULOUS FREQ - 
PRINT** 'IF YOU WISH TO USE THE EST. FREQ TYPE 1, ELSE TYPE 0 

READ* , I EST 

IF (I EST. EQ. 1) THEN 


SKIPC1 = 
SKIPC2 = 
SKIPC3 = 
DISCR1 = 
SKIPPR = 
PENDC1 = 
PENDC2 = 
PENDC3 = 
DISCR2 = 
PENDFR = 
IF 


, PENDFR 
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FIX = SKIPFR 
F1Y = SKIPFR 
F2X = PENDFR 
F2Y = PENDFR 
END IF 


C 


CALCULATE THE NOISELESS SIGNAL 
THE1X = TPIDT*F1X 
THE2X = TPIDT*F2X 
THE1Y = TPIDT*F1Y 
THE2Y = TPIDT*F2Y 
DO I = 1 , IPER 
II = 1-1 

X(I) = AOX + A1X*C0S(THE1X*I1+PHI1X) 
Y (I) = AOY + A1Y*C0S(THE1Y*I1+PHI1Y) 

END DO 


+A2X*COS (THE2X*I1+PHI2X) 
+A2Y*COS (THE2 Y*I1+PHI2Y) 


C LOOP TO INCORPORATE NOISE INTO 

IF (NOISE. EQ. 2) GO TO 1 
XRMS =0.0 
YRMS =0.0 
icount = 0.0 
XMAX = ABS (X ( 1) ) 

YMAX = ABS ( Y ( 1 ) ) 

DO I = 1 , IPER 

IF (ABS (X(I) ) .GT. XMAX) 

IF (ABS (Y (I) ) .GT .YMAX) 

2 vl = 2 . 0 * ranl(iseed) 

v2 = 2.0 * ranl(iseed) 
r = vl**2 + v2**2 
if (r.ge.l) go to 2 
fac = sqrt (-2 . 0*log(r) /r) 
xn = vl*fac*sd + x(i) 
yn = v 2 *f ac*sd + y(i) 

XRMS = XRMS + (XN - X(I))**2 
YRMS = YRMS + (YN - Y(I))**2 
X(I) = XN 
Y (I) = YN 
END DO 

XRMS = SQRT (XRMS/ (IPER+1) ) 

YRMS = SQRT (YRMS/ (IPER+1) ) 

XSNR = XMAX /XRMS 
YSNR = YMAX /YRMS 

PRINT*, 'XSNR = ' , XSNR, ' YSNR = ' , YSNR 
C LOOP TO SIMULATE DROPOUTS IN THE DATA 


THE DATA 


XMAX = ABS (X ( I) ) 
YMAX = ABS (Y (I) ) 
- 1.0 
- 1.0 


1 CONTINUE 

DROP1 = 1.0 - DROP 
DO I = 1 , IPER 

IF (ranl(ISEED) .GT.DROP1) THEN 



X ( I ) = o.o 
Y ( I ) = 0.0 
END IF 
END DO 

C WRITE OUT DATA TO IDFTXY.DAT 


DO I = 1 , IPER 

TIME = (1-1) *DT + TO 

WRITE ( 10, *) TIME, X (I) ,Y(I) , TLNGTH,MODE 
END DO 

CLOSE (10, STATUS = 'KEEP') 

END 


function ranl(idum) 

dimension r(97) . 

parameter (ml = 259200, ial = 7141, id 

parameter (m2 = 134456, ia2 = 8121, ic2 

parameter (m3 = 243000, ia3 = 4561, ic3 

data iff /0/ 

if (idum. It . 0 .or . if f . eq. 0) then 
iff = 1 


54773, rml = 1.0/ml) 
28411, rm2 = 1.0/m2) 
51349) 


ixl = mod ( icl — idum, ml) 

ixl = mod ( ial*ixl + icl, ml) 

ix2 = mod (ixl, m2) 

ixl = mod ( ial*ixl + icl, ml) 

ix3 = mod (ixl, m3) 

do j = 1,97 

ixl = mod(ial*ixl + icl, ml) 
ix2 = mod(ia2*ix2 + ic2,m2) 
r ( j ) = (float ( ixl) + f loat (ix2) *rm2) *rml 

end do 
idum = 1 
end if 


icl, ml) 
ic2 ,m2 ) 
ic3 ,m3) 


pause 


ixl = mod ( ial*ixl + 
ix2 = mod(ia2*ix2 + 
ix3 = mod ( ia3*ix3 + 
j = 1 + (97*ix3)/m3 
if ( j .gt. 97 .or. j . It. 1) 
rani = r ( j ) 

r ( j ) = (float ( ixl) + f loat (ix2) *rm2) *rml 


return 

end 
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C Program CRELIBR. FOR generates model data sets for the UNOMSC.FOR backup skip 
C rope observer program for the satellite at station 1 with a llbr atron compo- 
c nent in the tether. Five numbers are outputted per line of the data file 
c IDFTXY^DAT? time, x (roll) axis gyro rate, y (pitch) axis gyro rate, tether 
C length, and amcsmode. The user can opt to include noise and/or dropouts into 
C the^data . If desired, the program will estimate the skiprope and pendulous 
C frequencies as functions of the tether length. (Note: These estimated fre- 
C quencies are exactly that - estimates - and definitely should not be con- 

C strued as accurate values!) . , 

C The operator is prompted to input the following values: 
c time sampling interval, beginning time tag 

r number of data points, standard deviation, and random seed 

C (standard deviation is 10 X the nominal value of 2.8e-04) 

C (random seed should be a negative seed) 

C skiprope amplitude, frequency, and phase (for x axis) 

c skiprope amplitude, frequency, and phase (for y axis) 

c 1 or 2 to add noise or not, dropout percentage (0.0 to 1.0) 

r orb rate, pendulous amplitude, frequency, and phase (for x axis) 

c orb rate, pendulous amplitude, frequency, and phase (for y axis) 

C libration amplitude for x, libration amplitude for y, and x phase 

C (y ph ase is computed as x phase - 90.0) 

c tether length, amcsmode 

c 0 to use inputted values to generate model, 1 to estimate skiprope 

c and pendulous frequencies from the inputted tether length 

-C If noise is added to the data, the program prints to the operator the 

C values of the signal-to-noise ratios on each axis, XSNR and YSNR. 

c The signal-to-noise ratio is calculated as the maximum value divided 

C by the root mean square value for that axis. The noise is calculated 

' c using the Box-Muller method of generating Gaussian noise, and a portable 

C random number generator called rani found in "Numerical Recipes . 

C OPEN OUTPUT FILES ^ 'IDFTXY.DAT' (FOR UNOMSC.FOR) AND ' F_REC . DAT ' (FOR THE 
C ERROR CALCULATION PROGRAM VTBUSO.FOR) 

OPEN (9, FILE = 'F_REC.DAT' , STATUS = 'UNKNOWN') 

OPEN (10, FILE = 'IDFTXY.DAT' , STATUS = 'UNKNOWN') 

P PFAD INPUTS 

PRINT* , ' ENTER THE TIME SAMPLING INTERVAL DT AND BEGINNING TIME' 

DT TO 

PRINT* 'ENTER NUMBER OF TIME POINTS, STANDARD DEVIATION, RANDOM SEED' 
PRINT* (RANDOM SEED SHOULD BE A NEGATIVE NUMBER)' 

PRINT*, ' SKIPROPE AMPLITUDE ?, FREQUENCY ?, PHASE ? (FOR X)' 

PRINT* ,' SKIPROPE AMPLITUDE ?, FREQUENCY ?, PHASE ? (FOR Y) ' 

PRINT*, 'NOISE ? YES =1, NO = 2 ; % DROPOUT DESIRED (0.0 TO 1.0)' 

PRINT* ENTER ORB RATE, PENDULOUS AMP, FREQ, AND PHASE (FOR X)' 

PRINT* , ' ENTER ORB RATE, PENDULOUS AMP, FREQ, AND PHASE (FOR Y) ' 

READ* , A0Y , A2 Y , F2 Y , PHI2Y 
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PRINT* ENTER LIBRATION AMPLITUDE FOR X, Y, AND PHASE FOR X' 

~ READ* , ALX , ALY , PHILX 

PRINT* , ' ENTER THE TETHER LENGTH IN KILOMETERS AND AMCSMODE' 

READ* , TLNGTH , MODE 

C WRITE THE OUTPUT FILE 'F_REC.DAT' 

WRITE (9 , * ) NOISE, SD, ISEED 
_ WRITE (9 , *) A1X, FIX, PHI1X 

WRITE (9 , *) A1Y , F1Y , PHI1Y 

WRITE (9 , * ) A2X, F2X, PHI2X 

- C CONVERT PHASES FROM RADIANS TO DEGREES AND COMPUTE LIBRATION FREQUENCIES 

PHILY = PHILX - 90.0 
PI = 4 . 0*ATAN (1.0) 

CONVRT = PI/ 180.0 
PHI IX = PHI1X*C0NVRT 
PHI2X = PHI 2 X* CONVRT 
PHI1Y = PHI 1Y*C0NVRT 
PHI2Y = PHI2 Y*CONVRT 
PHILX = PH I LX* CONVRT 
PHILY = PHILY*CONVRT 

- TP IDT = 2 . 0*PI*DT 
FLX = 1/2713.0 
FLY = 1/3132.0 

C ESTIMATE THE SKIPROPE AND PENDULOUS FREQUENCIES FROM THE TETHER LENGTH 
SKIPC1 = 0.2119024 
SKIPC2 = -0.3571371 
“ SKIPC3 = 0.5309507 

DISCR1 = SKIPC2**2 - 4 . 0*SKIPC3 * (SKIPC1-TLNGTH) 

SKIPPR = 40 . 0* ( (-1) *SKIPC2 + SQRT (DISCR1) ) /SKIPC3 

- PENDC1 = 8 . 0952965E-02 
PENDC2 = 0.5571405 
PENDC3 = 0.2476189 

DISCR2 = PENDC2**2 - 4 . 0*PENDC3* (PENDC1-TLNGTH) 

- PENDFR = 0.01*(1.0+((-l)*PENDC2 + SQRT (DISCR2 ) ) / (2 . 0*PENDC3 ) ) 

IF (TLNGTH. LT. 1. 2) THEN 

PENDFR = PENDFR + 0.008 

- SKIPPR = 116 . 049*TLNGTH 
END IF 

SKIPFR = 1 0 /SKIPPR 

PRINT* ' EST. SKIPROPE FREQ = ', SKIPFR,' EST. PENDULOUS FREQ = ' , PENDFR 
PRINT* , ' IF YOU WISH TO USE THE EST. FREQ TYPE 1, ELSE TYPE 0' 

READ*, I EST 
IF (IEST.EQ.l) THEN 
“ FIX = SKIPFR 

F1Y = SKIPFR 
F2X = PENDFR 

- F2Y = PENDFR 
END IF 
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C CALCULATE THE NOISELESS SIGNAL 
THE1X = TPIDT*F1X 
THE2X = TPIDT*F2X 
THE1Y = TPIDT*F1Y 
THE2Y = TPIDT*F2Y 
THELX = TPIDT*FLX 
THELY = TPIDT*FLY 
DO I = 1 , IPER 

= aoX + A1X*C0S (THE1X*I1+PHI1X) +A2X*COS (THE2X*I1+PHI2X) 

g + ALX*COS(THELX*Il+PHILX) 

V y (I) = AOY + A1Y*C0S(THE1Y*I1+PHI1Y)+A2Y*C0S(THE2Y*I1+PHI2Y) 

& v + ALY*COS(THELY*Il+PHILY) 

END DO 

C LOOP TO INCORPORATE NOISE INTO THE DATA 

IF (NOISE. EQ. 2) GO TO 1 

XRMS =0.0 

YRMS =0.0 

icount = 0.0 

XMAX = ABS (X ( 1) ) 

YMAX = ABS ( Y ( 1) ) 

DO I = 1, IPER 

IF (ABS (X(I) ) .GT.XMAX) XMAX = ABS(X(I)) 

IF (ABS ( Y ( I) ) • GT . YMAX) YMAX = ABS(Y(I)) 

2 vl = 2.0 * ranl(iseed) - 1.0 

v2 = 2.0 * ranl(iseed) - 1.0 
r = vl**2 + v2**2 
if (r.ge. 1) go to 2 
fac = sqrt (-2 . 0*log(r) /r) 
xn = vl*fac*sd + x(i) 
yn = v 2 *fac*sd + y(i) 

XRMS = XRMS + (XN - X(I))**2 

YRMS = YRMS + (YN - Y(I))**2 

X(I) = XN 

Y (I) = YN 

END DO 

XRMS = SQRT (XRMS/ (IPER+1) ) 

YRMS = SQRT (YRMS/ (IPER+1) ) 

XSNR = XMAX /XRMS 
YSNR = YMAX /YRMS 

PRINT* , ' XSNR = ' t XSNR , ' YSNR = / , YSNR 
C LOOP TO SIMULATE DROPOUTS IN THE DATA 

1 CONTINUE 

DROP1 =1.0- DROP 
DO I = 1 , IPER 

IF (ranl(ISEED) .GT.DROP1) THEN 
X(I) = 0.0 
Y (I) = 0.0 
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END IF 
END DO 

C WRITE OUT DATA TO IDFTXY.DAT 

DO I = 1 , IPER 

TIME = ( 1-1) *DT + TO 

WRITE ( 10,*) TIME, X (I) , Y (I ) , TLNGTH , MODE 
END DO 

CLOSE (10, STATUS = 'KEEP' ) 

END 

function ranl(idum) 

dimension r(97) , , _ 

parameter (ml = 259200, ial — 7141, id — 54773, rml — 

parameter (m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 — 

parameter (m3 = 243000, ia3 = 4561, ic3 — 51349) 

data iff /0 / 

if (idum.lt.O.or.iff .eq.0) then 
iff = 1 

ixl = mod(icl - idum,ml) 

ixl = mod(ial*ixl + icl,ml) 

ix2 = mod (ixl, m2) 

ixl = mod(ial*ixl + icl,ml) 

ix3 = mod (ixl, m3) 

do j = 1,97 

ixl = mod ( ial*ixl + icl,ml) 
ix2 = mod(ia2*ix2 + ic2,m2) 
r ( j) = (float ( ixl) + f loat (ix2) *rm2) *rml 

end do 
idum = 1 
end if 

ixl = mod(ial*ixl + icl,ml) 
ix2 = mod ( ia2*ix2 + ic2,m2) 
ix3 = mod(ia3*ix3 + ic3,m3) 
j = 1 + (97*ix3)/m3 
if (j .gt.97.or. j . It. 1) pause 
rani = r(j) 

r ( j) = (float (ixl) ' + float (ix2) *rm2) *rml 

return 

end 


1 . 0/ml) 
1. 0/m2) 


9 



APPENDIX 2.B 


Results of ECR Verification Table 2. 


Filter Test Cases. 



This is the F_YAWMAN.DAT file for case 1 of the verification table 


midpoint time for this data window 

TIME TAG ON LAST POINT IN BUFFER IS 
REV LABEL POLARITY 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9439 
1.9439 
1.9439 
1.9439 
1.9439 


IS 

« 

# 


OF 


. 303104E+03 
. 102298E+04 
ROTATIONS START TIME 


.8 

.8 

.8 

.8 

.8 


1143.32 

1328.52 

1513.72 

1698.91 

1884.11 


This is the T_REPORT.DAT file for case 


1 of the verification table. 



INPUT 

OUTPUT 

ERROR 

Amp x 
Amp y 

.02000 

.02000 

.02001 

.01999 

-.041 % 
.063 % 

Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

. 000001 
.000000 

Phase x 
Phase y 

-107.0 

163.0 

-106.7 

163.3 

.32 

.35 


PHASE GRAD 
degs/cycle 


.054 

.010 


PHASE 6 T=15 MIN 
degs 


.58 

.39 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 

No noise for this run. . 

Stop index, start index, & number of point 


.000 .03125 10.0 

7.9 7.9 

= 593 1 593 


The following eleven pages represent the results of case 2 
of the verification table. The first ten pages give the results 
of the ten individual noise runs, with the F_YAWMAN.DAT file 
listed first and then the T_REPORT.DAT file. The last page lists 
the averaged results of the T_REPORT.DAT files. 



MIDPOINT TIME FOR THIS 
TIME TAG ON LAST POINT 
REV LABEL POLARITY 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


DATA WINDOW IS 
IN BUFFER IS 
RATE ( D/S ) 

1.9436 .8 

1.9436 *8 

1.9436 -8 

1.9436 -8 

1.9436 • 8 


. 303104E+03 
. 102298E+04 

# OF ROTATIONS START TIME 


1144.53 

1329.76 

1514.98 

1700.21 

1885.44 



INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

.02000 

.02000 

.01989 

.01955 

.565 % 
2.226 % 


Freq x 
Freq y 

.00540 

.00540 

.00540 

.00539 

.000005 Hz 
.000007 Hz 


Phase x 
Phase y 

-107.0 

163.0 

-107.5 

165.2 

.51 
2 . 16 

.325 

.489 


PHASE 8 T=15 MIN 
degs 


2.09 

4.53 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 

iseed = “1000 

Noise data: 1 sigma value (deg/sec) 

Stop index, start index, & number of point 


.000 

7.7 


.03125 

7.8 


10.0 


= 593 


000280 
1 593 


-midpoint time for this data window is 

SS 0„ LAST POINT IN BUFFER IS . 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.9452 
1.9452 
1.9452 
1.9452 
1.9452 


# OF ROTATIONS 
.8 
.8 
.8 
.8 
.8 


. 303104E+03 
. 102298E+04 


START TIME 
1143.61 
1328.68 
1513.75 
1698.82 
1883.89 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

.02000 

.02000 

.02005 

.01975 

-.230 % 
1.239 % 


.00540 

.00540 

.00540 

.00540 

.000003 Hz 
.000004 Hz 


-107.0 

163.0 

-107.4 

163.3 

.38 
. 34 

.169 

.279 


PHASE § T=15 MIN 
degs 


1.20 

1.70 


Pendulous data, amp, freq, phase _* 
Equivalent U & V deflections (meters) -. 
i coed = -2000 

Noise data: 1 sigma value (deg/sec): 
Stop index, start index, & number of point 


.000 

7.8 


.03125 

7.9 


10.0 


= 593 


.000280 
1 593 


MIDPOINT time for this data window is 
TIME^TAG ON LAST POINT IN BUFFER IS « 

n»mr / n / C \ 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S ) 
1.9444 
1.9444 
1.9444 
1.9444 
1.9444 


# OF 


. 303104E+03 
. 102298E+04 
ROTATIONS START TIME 


.8 

.8 

.8 

.8 

.8 


1143.05 

1328.20 

1513.35 

1698.50 

1883.65 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

.02000 

.02000 

.01974 

.02024 

1.322 % 
-1.204 % 


.00540 

.00540 

.00540 

.00540 

.000001 Hz 
.000001 Hz 

t 

-107.0 

163.0 

-106.9 

163.1 

. 10 
.08 

.084 

.056 


PHASE § T=15 MIN 
degs 


.50 

.35 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 

iseed = -3000 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point 593 


.000 

8.0 


.03125 

7.8 


000280 
1 593 


10.0 


midpoint time for this data window is 

SSTSS ON LAST POINT IN BUFFER IS . 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9440 
1 . 9440 
1.9440 
1.9440 
1 . 9440 


# OF ROTATIONS 
.8 
.8 
.8 
.8 
.8 


. 303 104E+03 
. 102298E+04 


START TIME 
1144.27 
1329.46 
1514.64 
1699.82 
1885.01 


INPUT 


.02000 

.02000 

.00540 

.00540 

-107.0 

163.0 


OUTPUT 


.02024 

.01972 

.00540 

.00540 

-107.3 

164.2 


ERROR 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 

Pendulous data, amp, freq, phase =: 

Equivalent U & V deflections (meters) -. 
iceed = -4000 j . . 

Noise data: 1 sigma value (deg/sec): 

Stop index, start index, 6 number o£ point 


•1.224 % 
1.384 % 

.000002 Hz 
.000002 Hz 

.33 

1.20 


PHASE GRAD 
degs/ cycle 


, 139 
, 132 


PHASE 8 T=15 MIN 
degs 


1.01 

1.84 


.000 

7.8 


.03125 

8.0 


10.0 


= 593 


.000280 
1 593 



midpoint time for this data window is 

TIME TAG ON LAST POINT IN BUFFER IS • 

* * ■» mn / n / P \ 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE ( D/S ) 
1.9405 
1.9405 
1.9405 
1.9405 
1.9405 


. 303104E+03 
. 102298E+04 
# OF ROTATIONS START TIME 
.8 1145.21 

.8 1330.73 

.8 1516.24 

.8 1701.76 

.8 1887.28 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

.02000 

.02000 

.01991 

.02040 

.00540 

.00540 

.00539 

.00539 

-107.0 

163.0 

-107.5 

164.8 


ERROR 


.448 % 
-2.024 % 

.000009 Hz 
.000010 Hz 

.52 

1.85 


Pendulous data, amp, freq, phase _* 

Equivalent U & V deflections (meters) 
is6@d = -5000 

Noise data: 1 sigma value (deg/sec): 

Stop index, start index, & number of point 


PHASE GRAD 
degs/cycle 


.610 

.679 


PHASE 0 T=15 MIN 
degs 


3.48 

5.15 


.000 

8.0 


.03125 

7.8 


10.0 


= 593 


,000280 
1 593 



MIDPOINT TIME FOR THIS DATA WINDOW IS 
SSfwS ON LAST POINT IN BUFFER IS = 


REV LABEL 
0 
1 
2 

3 

4 


Amp x 
Amp y 


Freq x 
Freq y 


Phase x 
Phase y 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


INPUT 


.02000 

.02000 


.00540 

.00540 


-107.0 

163.0 


RATE (D/S) 
1.9412 
1.9412 
1.9412 
1.9412 
1.9412 


OUTPUT 


.02028 

.02004 


.00540 

.00538 


-106 . 3 
164.8 


: . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 


.8 

.8 

.8 

.8 

.8 


ERROR 


■1.423 % 
-.193 % 


.73 

1.80 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 
tcood = -6000 

Noise data: 1 sigma value (deg/sec): 

Stop index, start index, & number of point 


1144.85 

1330.31 

1515.76 

1701.22 

1886.68 


PHASE GRAD 
degs/cycle 


PHASE @ T=15 MIb 
degs 


.000001 Hz 
.000016 Hz 


.044 

1.097 


.95 
7 . 13 


.000 

7.9 


.03125 10.0 

8.0 


= 593 


,000280 
1 593 



MIDPOINT time for this data window is 
ON LAST POINT IN BUFFER IS • 

1 » mn / n / C 1 \ 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9405 
1.9405 
1.9405 
1.9405 
1.9405 


. 303104E+03 
. 102298E+04 
# OF ROTATIONS START TIME 
.8 1144.19 

8 1329.71 

.8 1515.23 

.8 1700.75 

8 1886.27 


INPUT 

OUTPUT 

ERROR 

.02000 

.01992 

. 383 % 

.02000 

.01996 

. 184 % 

. 00540 

.00539 

.000008 Hz 

.00540 

.00539 

.000012 Hz 

-107.0 

-105.6 

1.42 

163.0 

164.8 

1 .79 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 

Pendulous data, amp, freq, phase =: 

Equivalent U & V deflections (meters) 
iseed = -7000 

Noise data: 1 sigma value (deg/sec): 

Stop index, start index, & number of point 


PHASE GRAD 
degs /cycle 


.502 

.791 


PHASE 0 T=15 MIN 
degs 


3.86 

5.64 


.000 

7.9 


.03125 

7.8 


10.0 


= 593 


.000280 
1 593 


midpoint time for this data window is 

TIME TAG ON LAST POINT IN BUFFER IS • 
REV LABEL POLARITY 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE ( D/S ) 
1.9454 
1.9454 
1.9454 
1.9454 
1.9454 


. 303104E+03 


. 102298E+04 


ROTATIONS 

START TIME 

.8 

1142.47 

.8 

1327.52 

.8 

1512.57 

.8 

1697.62 

.8 

1882 .67 



INPUT 

OUTPUT 

ERROR 

Amp x 
Amp y 

.02000 

.02000 

.01985 

.02034 

.737 % 
-1.681 % 

Freq x 
Freq y 

.00540 

.00540 

.00539 

.00541 

.000005 Hz 
.000013 Hz 

Phase x 
Phase y 

-107.0 

163.0 

-107.4 

162.0 

.39 

1.01 


PHASE GRAD 
degs/cycle 


.357 

.880 


PHASE § T=15 MIN 
degs 


2.12 

5.28 


Pendulous data, amp, freq, phase _• 

Equivalent U & V deflections (meters) -. 

{ seed = -8000 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point - 593 


.000 

8.0 


,03125 

.8 


10.0 


.000280 
1 593 



MIDPOINT TIME FOR THIS 
TIME TAG ON LAST POINT 
REV LABEL POLARITY 

0 POSITIVE 

1 POSITIVE 

2 POSITIVE 

3 POSITIVE 

4 POSITIVE 

DATA WINDOW 
IN BUFFER IS 
RATE (D/S) 
1.9464 
1.9464 
1.9464 
1.9464 
1.9464 

IS : . 303104E+03 

: . 102298E+04 

# OF ROTATIONS START TIME 
.8 1142.96 

.8 1327.91 

.8 1512.87 

.8 1697.83 

.8 1882.79 



INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

PHASE § T=15 
degs 

Amp x 
Amp y 

.02000 

.02000 

.01977 

.02004 

1.143 % 
-.209 % 



Freq x 
Freq y 

.00540 

.00540 

.00541 

.00541 

.000008 Hz 
.000006 Hz 



Phase x 
Phase y 

-107.0 

163.0 

-108.3 

163.2 

1.29 

.20 

.501 

.379 

3.73 

2.04 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 
iseed — -9000 

Noise data: 1 sigma value (deg/sec) 

Stop index, start index, & number of point 


.000 

7.9 


.03125 

7.8 


10.0 


= 593 


.000280 
1 593 



MTDPOINT time for this data window is 
TAG ON LAST POINT IN BUFFER IS : 
REV LABEL POLARITY 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9456 
1.9456 
1.9456 
1.9456 
1.9456 


# OF 


. 303104E+03 
. 102298E+04 
ROTATIONS START TIME 
8 1142.38 

.8 1327.41 

.8 1512.45 

.8 1697.48 

o 1882.52 



INPUT 

OUTPUT 

ERROR 

Amp x 
Amp y 

.02000 

.02000 

.01999 

.02025 

.069 % 
-1.240 % 

Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000004 Hz 
.000005 Hz 

Phase x 
Phase y 

-107.0 

163.0 

-106.5 

163.3 

.45 

.28 


PHASE GRAD 
degs/cycle 


.256 

.328 


PHASE § T=15 MIN 
degs 


1.70 

1.87 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) -. 

tseed = -10000 j , 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point - 593 


.000 

8.0 


.03125 

7.8 


10.0 


.000280 
1 593 




INPUT 

OUTPUT 

ERROR 

Amp x 
Amp y 

.02000 

.02000 

.01996 

.02003 

.179 % 
-.152 % 

Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000005 Hz 
.000008 Hz 

Phase x 
Phase y 

-107.0 

163.0 

-107.1 

163.9 

.61 

1.07 


Pendulous data/ amp/ freq, phase 
Fauivalent U & V deflections (meters) . 

q Noise data: 1 sigma value 

Stop index, start index, & number of pom 


PHASE GRAD 
degs/cycle 


PHASE @ T=15 MIN 
degs 


.299 2 * 06 

.511 3 - 55 

.000 .03125 10.0 

7.9 7.8 

.000280 
593 1 593 



i n f the verification table. 

This is the F_YWAMAN.DAT me for case 3 


rn n mtiT c DATA WINDOW IS 
MIDPOINT TIME FO TH BUFFER is : 

TI ME TAG ON LAST POINT /g) # 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.6581 
1.6581 
1.6581 
1.6581 
1.6581 


. 303104E+03 
. 102298E+04 

OF ROTATIONS START TIME 

.9 
.9 
.9 
.9 
.9 


1028.70 

1245.81 

1462.92 

1680.03 

1897.14 


.cm ■Fr'.r- rase 3 of the verification table 
This is the T_REPORT.DAT file for case 


Amp x 
Amp y 

Freq x 
Freq y 

phase x 
Phase y 


INPUT 


.02000 

.02000 

.00460 

.00460 

50.0 

-40.0 


OUTPUT 


.02002 

.01995 

.00460 

.00461 

49.2 

-41.3 


ERROR 


-.099 % 
.269 % 

,000003 Hz 
.000009 Hz 

.80 

1.28 


PHASE GRAD 
degs/cycle 


PHASE S T=15 MIN 
degs 


.226 

.705 


1.74 

4.20 


Equivalent'u^ ’ V*S i lefuinf titers , "= I 

Stop'index^ start Index , . nu^her of point - 593 


.500 

9.2 


.03125 

9.2 


60.0 


1 593 



The following eleven pages represent the results of case 4 
of the"verif ication table." The ten P ^ the results 

i-hP ten individual noise runs, with the F_yawmaw.uai 
Usted fi£s£ and then the T_REPORT.DAT file. The last page lists 
the averaged results of the T_REPORT.DAT files. 



_ T1U(C , voR THIS DATA WINDOW IS 

as " ; 

nrvr T rn T \7V 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


1.6565 

1.6565 

1.6565 

1.6565 

1.6565 


. 303104E+03 
. 102298E+04 


rotations 

START TIME 

# 9 

1029.51 

. 9 

1246.83 

. 9 

1464.15 

, 9 

1681.47 

.9 

1898.79 


input output 


ERROR 


PHASE GRAD PHASE § T-15 MIN 

degs/cycle degs 


Amp x 
Amp y 

Freq x 
Freq y 


.02000 

.02000 

.00460 

.00460 


.01989 

.02014 

.00460 

.00460 


Phase x 
Phase y 


50.0 

-40.0 


49.4 

-40.7 


.556 % 
-.676 % 

.000001 Hz 
.000004 Hz 

.61 

.72 


Pendulous data, amp, freq, . \ 

Equivalent U i V deflections (meters) 

ise ® d ~ data- 0 1 sigma value (deg/sec): 
StopTnde*fs«rt indel „ number of point - 


.044 

.281 

.500 .03125 

9.3 9.2 


.000280 
593 1 593 


.79 

1.88 


60.0 



midpoint time FOR THIS 
time tag on last point 

REV LABEL POLARITY 
0 


1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


DATA WINDOW IS 
IN BUFFER IS : 
RATE(D/S) # 

1.6618 
1.6618 
1.6618 
1.6618 
1.6618 


. 303104E+03 
102298E+04 


OF ROTATIONS 
.9 
.9 
.9 
.9 
.9 


START TIME 
1026.85 
1243.48 
1460.11 
1676.75 
1893.38 


INPUT OUTPUT 


ERROR 


PHASE GRAD PHASE G T-15 MIN 

degs /cycle degs 


Amp x 
Amp y 


.02000 

.02000 


Freq x 
Freq y 


.00460 

.00460 


Phase x 3U •" 

Phase y -40.0 


.02009 * 

.02020 -1.014 % 


.00462 

.00462 


.000016 Hz 
.000016 Hz 


48.1 i- 94 

-42.5 2.49 


Pendulous data, amp, freq, P h ^® . J 

Equivalent U & V deflections (meters) • 

lse f " dater 0000 1 sigma value (deg/sec) s 
S top° indexes tart indel a nu.ber o£ point - 


1.244 

1.276 

.500 .03125 

9.3 9.2 


.000280 
593 1 593 


7.09 

7.77 


60.0 



MTDPOINT time for this data window is 

TIME°tAG ON LAST POINT IK* IS = 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.6579 
1.6579 
1.6579 
1.6579 
1.6579 


.303104E+03 
. 102298E+04 
# OF ROTATIONS START TIME 
9 1029.86 

' 9 1247.00 

* 9 1464.15 

* 9 1681.30 

*g 1898.44 


input 

OUTPUT 

ERROR 

.02000 

.02000 

.01980 

.01971 

.979 % 
1.441 % 

.00460 

.00460 

.00460 

.00461 

.000000 Hz 
.000010 Hz 

50.0 

-40.0 

48.8 

-40.5 

1.18 

.51 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 

Pendulous data, amp, freq, P^e . 
Equivalent U 6 V deflections (meters) . 

iseed - ^ sigma value (deg/sec) '• 

Stop°ln*e* =tart indeS, 4 numter of point 


PHASE GRAD 
degs /cycle 


PHASE 0 T=15 MIN 
degs 


.001 

.814 


1.19 

3.88 


.500 
l . 1 


.03125 
9 . 1 


60.0 


= 593 


.000280 
1 593 



mTMP for THIS DATA WINDOW IS 
midpoint time FOR THi BUFF er IS : 

TIME TAG ON LAST POINT /g) # 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE ( D/S ) 
1.6587 
1.6587 
1.6587 
1.6587 
1.6587 


. 303104E+03 
. 102298E+04 

OF ROTATIONS START TIME 


.9 

.9 

.9 

.9 

.9 


1029.43 

1246.47 

1463.51 

1680.54 

1897.58 


input OUTPUT 


ERROR 


PHASE GRAD 
degs/cycle 


PHASE § T=15 MIN 
degs 


Amp x 
Amp y 


.02000 

.02000 


Freq x 
Freq y 


. 00460 
.00460 


Phase x 
phase y 


50.0 

-40.0 


.01973 

.02016 

.00461 

.00460 

47.7 

-40.8 


1.347 % 
-.799 % 

.000010 Hz 
.000005 Hz 

2.29 

.80 


iseed = -400000 valu e (deg / s ec) : 

HO -ndet Siart index, & number of point - 
Stop index, swll 


.819 

.357 

.500 .03125 

9.3 9 • 1 

.000280 
593 1 593 


5.68 

2.28 


60.0 



-SF sHiir” 

REV LABEL POLARITY 

0 POSITIVE 

1 POSITIVE 

2 POSITIVE 

3 POSITIVE 

4 POSITIVE 


. 303104E+03 
. 102298E+04 


RATE (D/S) 
1.6575 
1.6575 
1.6575 
1.6575 
1.6575 


ROTATIONS 

START TIME 

. 9 

1029.04 

. 9 

1246.24 

1 9 

1463.44 

. 9 

1680.63 

.9 

1897.83 


INPUT 

OUTPUT 

ERROR 

.02000 

.02000 

.02001 

.02022 

-.071 % 
-1.082 % 

.00460 

.00460 

.00461 

.00460 

.000007 

.000001 

50.0 

-40.0 

49.0 

-40.8 

1.04 

.83 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 

Pendulous data, amp, freq, P^se * 
Equivalent U & V deflections (meters) . 

iSe S d |^o data^° 00 °l Sigma value (deg/sec) : 
s top° indexes tart indel . number of point 


PHASE GRAD 
degs/cycle 


PHASE @ T=15 MIN 
degs 


.570 

.068 


3.40 

1.11 


.500 

9.3 


.03125 

9.2 


60.0 


= 593 


.000280 
1 593 



midpoint time for this data window !S 

time TAG ON LAST POINT IN BUFFER . 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.6582 
1.6582 
1.6582 
1.6582 
1.6582 


. 303104E+03 
. 102298E+04 

OF ROTATIONS START TIME 


.9 

.9 

.9 

.9 

.9 


1028.65 

1245.75 

1462.85 

1679.95 

1897.05 


INPUT OUTPUT 


ERROR PHASE GRAD 

degs/cycle 


PHASE 8 T=15 MIN 
degs 


Amp x 
Amp y 

Freq x 
Freq y 


.02000 

.02000 

.00460 

.00460 


.02017 

.01987 

.00461 

.00460 


Phase x 50 
Phase y -40.0 


49.5 

-40.1 


-.845 % 
.649 % 

.000013 Hz 
.000000 Hz 


.49 

.08 


Pendulous data, amp, fr ®<3' ph ^®,. prs \ = ’. 
Equivalent U & V deflections (meters) . 

lse ? d ( ! data^ 000<> l sigma value (deg/sec) : 
Stop° indexes tart US. « number of point - 


.996 

.027 

.500 .03125 

9.2 9.3 

.000280 
593 1 593 


4.62 

.20 


60.0 



midpoint time for this data window is 

S3TSS ON LAST POINT IN BUFFER IS = 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.6577 
1.6577 
1.6577 
1.6577 
1.6577 


. 303104E+03 
. 102298E+04 
# OF ROTATIONS START TIME 
g 1029.96 

* 9 1247.13 

'9 1464.30 

'9 1681.47 

* 9 1898.65 


INPUT 

OUTPUT 

error 

.02000 

.02000 

.01965 

.02015 

1.757 % 
-.735 % 

.00460 

.00460 

.00461 

.00460 

.000012 

.000003 

50.0 

-40.0 

47.9 

-40.3 

2.06 

.26 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 

Pendulous data, amp, freq, phase 
Equivalent U s V deflections (meters) 

Iseed * :'! 00000 1 slgma value (deg/sec) : 

Stop° indexes tart inCel & number of point 


PHASE GRAD 
degs/cycle 


PHASE @ T= 15 MIN 
degs 


.930 

.203 


5.91 

1.11 


.500 

9.3 


.03125 60.0 

9 . 1 


593 


.000280 
1 593 



™“ SSWS! SMT-” 

REV LABEL POLAR! 1* 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.6561 
1.6561 
1.6561 
1.6561 
1.6561 


. 303104E+03 


a x \J *< *» - r 

ROTATIONS 

START TIME 

g 

1029.72 

9 

1247.10 

9 

1464.47 

9 

1681.85 

.9 

1899.22 


input 

OUTPUT 

ERROR 

.02000 

.02000 

.01974 

.01975 

1.309 % 
1.247 % 

.00460 

.00460 

.00460 

.00460 

.000000 Hz 
.000001 Hz 

50.0 

-40.0 

49.6 

-41.4 

.39 

1.42 


Amp x 
Amp y 

Freq x 
Freq Y 

phase x 
Phase y 

Pendulous data, amp, freq, P^se . 

Equivalent U s V deflectrons (meters) . 

lse » d -ne data® 00 1 sigma value (deg/sec) : 
Stop 0 indexes tart indel a number of pornt 


PHASE GRAD 
degs/cycle 


PHASE § T=15 MIN 
degs 


.030 

.083 


.51 

1.76 


.500 

9.1 


.03125 60.0 

9.1 


593 


.000280 
1 593 



MIDPOINT time FOR this oata WINDOW^S 

TIME TAG ON LAST POINT I" # 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.6571 
1.6571 
1.6571 
1.6571 
1.6571 


. 303104E+03 
. 102298E+04 

OF ROTATIONS START TIME 


.9 

.9 

.9 

.9 

.9 


1029.24 
1246.49 
1463.74 
1681.00 

1898.25 


INPUT OUTPUT 


ERROR 


PHASE GRAD PHASE 6 T=15 MIN 

degs/cycle degs 


Amp x 
Amp y 


.02000 

.02000 


Freq x 
Freq y 


.00460 

.00460 


Phase x 
Phase y 


50.0 

-40.0 


.02022 

.01993 


-1.115 % 
.348 % 


. 00460 
.00461 


.000001 Hz 
.000007 Hz 


49.8 

-40.8 


.22 

.84 


iS6 M d ^ e data- 0000 l sigma value (deg/sec) : 
Stoplndex 'tart indei S nu^er of point - 


.111 

.575 

.500 .03125 

9.2 9-3 

.000280 
593 1 593 


.68 

3.22 


60.0 



MIDPOINT TIME FO r this ^TgyppE^IS : 
time TAG on LAST POINT I* 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.6632 
1.6632 
1.6632 
1.6632 
1.6632 


. 303104E+03 
. 102298E+04 

# OF ROTATIONS START TIME 


.9 

.9 

.9 

.9 

.9 


1026 . 13 
1242.58 
1459.03 
1675.48 
1891.93 


INPUT OUTPUT 


ERROR 


PHASE GRAD 
degs/cycle 


PHASE 6 T=15 MIN 
degs 


Amp x 
Amp y 


.02000 

.02000 


Freq x 
Freq Y 


.00460 

.00460 


Phase x 
Phase y 


50.0 

-40.0 


.02002 

.02016 

.00463 

.00461 

47 .5 
-40.7 


-.113 % 
-.806 % 

.000026 Hz 
.000014 Hz 

2.53 

.66 


lSeed - = dati- 00 1 sigma ™lue ( de 9 /se< T> =. 
5 top° indexes tart Indel . nu^er of pornt - 


2.070 

1.062 

.500 .03125 

9.3 9.2 


.000280 
593 1 593 


11.10 

5.06 


60.0 



Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 


.02000 

.02000 

.00460 

.00460 

50.0 

-40.0 


OUTPUT 


.01993 

.02003 

.00461 

.00461 

48.7 

-40.9 


ERROR 


.336 % 
-.143 % 

. 000009 
. 000006 

1.28 

.86 


Hz 

Hz 


P " n ? U iCt d S“' 'vThec^inf titers) - 

Equivalent U sigma value (deg/sec) : 

Sto^index start indel, 4 nu^er of point 


PHASE GRAD PHASE @ T=15 Mil 

degs/cycle degs 


.682 

.475 


4.10 

2.83 


.500 .03125 60.0 

9.2 S.2 

.000280 
593 1 593 



The following eleven pages represent the results of case 5 
of the verification table. The first ten pages give the results 
of the ten individual noise runs, with the F_YAWMAN.DAT file 
listed first and then the T_REPORT.DAT file. The last page lists 
the averaged results of the T_REPORT.DAT files. 


MIDPOINT TIME FOR THIS DATA WINDOW IS 
TIME TAG ON LAST POINT IN BUFFER IS ■ 

X n»mr« / n / 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9436 
1.9436 
1.9436 
1.9436 
1.9436 


. 303104E + 03 
. 102298E+04 
# OF ROTATIONS START TIME 
5 9 1190.58 

5.9 1375.80 

5.9 1561.02 

5.9 1746.24 

5.9 1931.46 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

. 15000 

. 14993 

.048 % 

. 15000 

. 15022 

-.146 % 

.00540 

.00540 

.000001 

.00540 

.00540 

.000001 

160.0 

160.5 

.53 

70.0 

70.4 

. 37 


PHASE GRAD 
degs/cycle 


PHASE @ T=15 MI 
degs 


.073 

.068 


.89 

.70 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 

iseed = -11000 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point - 593 


.000 

59.0 


.03125 

58.9 


50.0 


.000280 
1 593 



midpoint time for this data window is 
S55°3S ON LAST POINT IN BUFFER IS . 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.9437 
1.9437 
1.9437 
1.9437 
1.9437 


: . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 
5 9 1190.56 

5.9 1375.78 

5.9 1560.99 

5.9 1746.21 

5 9 1931.43 



INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

. 15000 
. 15000 

. 15003 
. 15015 

-.020 % 
-.100 % 


Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000000 Hz 
.000002 Hz 


Phase x 
Phase y 

160.0 

70.0 

160.4 

70.5 

.42 

.53 

.026 

.102 


PHASE @ T=15 MIN 
degs 


.54 

1.03 


Pendulous data, amp, freq, phase _* 

Equivalent U & V deflections (meters) -. 

Iseed = -12000 j . 

Noise data: 1 sigma value (deg/sec) 

Stop index, start index, & number of point 593 


.000 

59.0 


.03125 

59.0 


50.0 


.000280 
1 593 



- MIDPOINT 

TIME FOR THIS 

DATA WINDOW 

TIME TAG 

ON LAST POINT 

IN BUFFER IS 

REV LABEL POLARITY 

RATE (D/S) 

0 

POSITIVE 

1.9434 

1 

POSITIVE 

1.9434 

2 

POSITIVE 

1.9434 

3 

POSITIVE 

1.9434 

4 

POSITIVE 

1.9434 


: , 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 


5.9 

5.9 

5.9 

5.9 

5.9 


1190.72 

1375.97 

1561.22 

1746.46 

1931.71 



INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs /cycle 

Amp x 
Amp y 

. 15000 
. 15000 

. 15011 
. 15012 

-.074 % 
-.082 % 


Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000001 Hz 
.000002 Hz 


Phase x 
Phase y 

160.0 

70.0 

160.4 

70.6 

.43 

.59 

.074 

.164 


PHASE 0 T=15 MIN 
degs 


.79 

1.39 


Pendulous data, amp, freq, phase • 
Equivalent U & V deflections (meters) — . 
iseed = “13000 

Noise data: 1 sigma value (deg/sec) : 
Stop index, start index, & number of point 


.000 

59.0 


.03125 

59.0 


50.0 


= 593 


.000280 
1 593 



midpoint TIME for THIS °£TA WINDOWS 

TIME TAG ON last POINT X 
REV LABEL 
0 
1 
2 

3 

4 


. 303104E+03 


i 

POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.9436 
1.9436 
1.9436 
1.9436 
1.9436 


, x V *• ** ^ 

ROTATIONS 
c, 9 

START TIME 
1190.57 

S 9 

1375.79 

5.9 

c Q 

1561.02 

1746.24 

5.9 

1931.46 


INPUT OUTPUT 


ERROR 


PHASE GRAD 
degs/cycle 


PHASE 0 T=15 MIN 
degs 


Amp x 
Amp y 


. 15000 
. 15000 


Freq x 
Freq Y 


.00540 

.00540 


Phase x 16 2‘? 

Phase y 70 -° 


. 14979 
. 14998 

.00540 

.00540 

160.5 

70.2 


.140 % 
.013 % 

.000001 Hz 
.000001 Hz 

.47 

.23 


Equivalent^U^ ' t-ters , ‘-i 
iseed " j_ ta T 1400 i sigma value (deg/sec): 
Stop° indexes tart Ul, A number of point - 


.093 

.046 

.000 .03125 

59.0 58.9 

.000280 
593 1 593 


.92 

.46 


50.0 



- MIDPOINT TIME FOR THIS 
TIME TAG ON LAST POINT 
REV LABEL POLARITY 

0 POSITIVE 

1 POSITIVE 

2 POSITIVE 

3 POSITIVE 

— 4 POSITIVE 


DATA WINDOW 

IS 

• 

• 

. 303104E+03 

IN BUFFER IS 

• 

• 


. 102298E+04 

RATE ( D/S ) 

# 

OF ROTATIONS START TIME 

1.9438 


5.9 

1190.49 

1.9438 


5.9 

1375.69 

1.9438 


5.9 

1560.90 

1.9438 


5.9 

1746.10 

1.9438 


5.9 

1931.31 



INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

.15000 
. 15000 

. 15010 
. 15037 

-.068 % 
-.246 % 


Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000000 Hz 
.000001 Hz 


Phase x 
Phase y 

160.0 

70.0 

160.2 

70.6 

.24 

.59 

.013 

.094 


PHASE @ T=15 MIN 
degs 


.31 

1.05 


Pendulous data, amp, freq, phase _* 

Equivalent U & V deflections (meters) 
iseed = -15000 

- Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point 


.000 .03125 50.0 

59.1 59.0 

.000280 
593 1 593 




MIDPOINT TIME FOR THIS DATA WINDOW IS 

~ TIME tag on last point in buffer is • 

REV LABEL POLARITY 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.9446 
1.9446 
1.9446 
1.9446 
1.9446 


: . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 


5.9 

5.9 

5.9 

5.9 

5.9 


1190.04 

1375.18 

1560.31 

1745.44 

1930.57 



INPUT 

OUTPUT 

ERROR 


PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

. 15000 
.15000 

. 14995 
. 15009 

.031 % 
-.060 % 



Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000003 

.000001 

Hz 

Hz 


Phase x 
Phase y 

160.0 

70.0 

160.0 

70.1 

.01 

.13 


. 168 
.043 


PHASE g T=15 MIN 
degs 


.82 

.34 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 
iseed = “16000 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point 


.000 

59.0 


.03125 

58.9 


50.0 


= 593 


000280 
1 593 



midpoint time for this data window is 

TIME TAG ON LAST POINT IN BUFFER IS * 

* * * *“ _ _ n » mD / n / O \ 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9443 
1.9443 
1.9443 
1.9443 
1.9443 


# OF ROTATIONS 
5.9 
5.9 
5.9 
5.9 
5.9 


. 303104E+03 
. 102298E+04 


START TIME 
1190.21 
1375.37 
1560.53 
1745.69 
1930.84 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

. 15000 
. 15000 

. 14987 
. 14994 

.084 % 
.040 % 


.00540 

.00540 

.00540 

.00540 

.000000 Hz 
.000001 Hz 


160.0 

70.0 

160.5 

70.1 

.50 
. 14 

.012 

.090 


PHASE § T=15 MIN 
degs 


.56 

.58 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 

.Led = -17000 

Noise data: 1 sigma value (deg/sec): 

Stop index, start index, & number of point 


.000 .03125 50.0 

58.9 58.9 


= 593 


.000280 
1 593 



midpoint time for this data window is 
S£°S5 ON LAST POINT IN BUFFER IS = 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9437 
1.9437 
1.9437 
1.9437 
1.9437 


. . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 
5 9 1190.51 

5 *9 1375.72 

5 ' 9 1560.93 

5 * 9 1746.14 

5*9 1931.35 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

. 15000 
. 15000 

. 14977 
.14986 

.152 % 
.096 % 

.00540 

.00540 

.00540 

.00540 

.000000 

.000001 

160.0 

70.0 

160.5 

70.2 

.55 
. 19 


PHASE GRAD 
degs/cycle 


PHASE @ T=15 MIN 
degs 


.029 

.067 


.69 

.51 


Pendulous data/ amp/ freq/ phase 
Equivalent U & V deflections (meters) . 

iseed = ” 18 °°i siama value (deg/sec) : 

Stop° indexes tart indel * number of point - 593 


.000 

58.9 


.03125 

58.9 


50.0 


.000280 
1 593 



midpoint time for this data window is 
time tag on last point in buffer is • 

REV LABEL POLARITY — 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9440 
1.9440 
1.9440 
1.9440 
1.9440 


: . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 
5.9 1190.34 

5.9 1375.52 

5.9 1560.70 

5.9 1745.88 

5.9 1931.06 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

. 15000 
. 15000 

. 14949 
. 14994 

.343 % 
.038 % 


.00540 

.00540 

.00540 

.00540 

.000001 Hz 
.000001 Hz 


160.0 

70.0 

160.2 

70.3 

.22 

.27 

.067 

.049 


PHASE 6 T=15 MIN 
degs 


.55 

.51 


Pendulous data, amp, freq, phase _• 

Equivalent U & V deflections (meters) 

1 seed = -19000 

Noise data: 1 sigma value (deg/sec): 

Stop index, start index, & number of point 


.000 .03125 50.0 

58.9 58.7 


= 593 


,000280 
1 593 



midpoint time for this data window is 
time°tag on last point in BUFFER is = 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

positive 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE(D/S) 
1.9439 
1.9439 
1.9439 
1.9439 
1.9439 


: . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 
5 9 1190.42 

5 ' 9 1375.61 

5 * 9 1560.80 

5 * 9 1746.00 

5 ’ 9 1931.19 



INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

. 15000 
. 15000 

. 14975 
. 15013 

.168 % 
-.089 % 


Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000001 Hz 
.000001 Hz 


Phase x 
Phase y 

160.0 

70.0 

160.2 

70.2 

. 19 
.25 

.045 

.080 


PHASE @ T=15 MIN 
degs 


.41 

.64 


Pendulous data/ amp, frog/ phase 
Equivalent U & V deflections (meters) -. 

1 seed = -20000 . . . 

Noise data: 1 sigma value ( de 9 / s® < ?) • 

Stop index, start index, & number of point - 593 


.000 

59.0 


.03125 

58.8 


50.0 


.000280 
1 593 




INPUT 

OUTPUT 

ERROR 

Amp x 
Amp y 

. 15000 
. 15000 

. 14988 
. 15008 

.080 % 
-.054 % 

Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000001 

.000001 

Phase x 
Phase y 

160.0 

70.0 

160.3 

70.3 

.36 

.33 


PHASE GRAD 
degs/cycle 


.060 

.080 


PHASE 8 T=15 MIN 
degs 


.65 

.72 


Pendulous data, amp, freq, phase : 

Equivalent U & V deflections (meters) 

Noise data: 1 sigma value (deg/sec). 

Stop index, start index, & number of point 


.000 

59.0 

= 593 


.03125 
58.9 
000280 
1 593 


50.0 



The following eleven pages represent the results of case 6 
of the verification table. The first ten pages give the results 
of the ten individual noise runs, with the F_YAWMAN.DAT file 
listed first and then the T_REPORT.DAT file. The last page lists 
the averaged results of the T REPORT.DAT files. 



MIDPOINT TIME FOR THIS 
TIME TAG ON LAST POINT 
REV LABEL POLARITY 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


DATA WINDOW IS 
IN BUFFER IS : 
RATE (D/S) # 

1.9441 
1.9441 
1.9441 
1.9441 
1.9441 


. 303104E+03 
. 102298E+04 


OF ROTATIONS 
5.9 
5.9 
5.9 
5.9 
5.9 


START TIME 
1190.31 
1375.49 
1560.67 
1745.84 
1931.02 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

. 15000 
. 15000 

. 14984 
. 15014 

.104 % 
-.096 % 


.00540 

.00540 

.00540 

.00540 

.000000 Hz 
.000001 Hz 


160.0 

70.0 

160.4 

70.1 

.45 

.05 

.014 

.048 


PHASE @ T=15 MIN 
degs 


.51 

.28 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) -. 
iseed = -110000 

Noise data: 1 sigma value (deg/sec): 

Stop index, start index, & number of point 


.500 

59.0 


.03125 

58.9 


50.0 


= 593 


,000280 
1 593 



ssrs; S-E.MS? Sir “ 

REV LABEL POLARITY RATE (D/S) 


0 

1 

2 

3 

4 


POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9436 
1.9436 
1.9436 
1.9436 
1.9436 


. 303104E+03 
. 102298E+04 

# OF ROTATIONS START TIME 

* % q 1190.61 

eg 1375.84 

eg 1561.07 

e'g 1746.29 

5 # g 1931.52 


INPUT 


OUTPUT 


ERROR 


PHASE GRAD 
degs/cycle 


PHASE § T=15 MI 
degs 


Amp x 
Amp y 

Freq x 
Freq y 

phase x 
Phase y 


, 15000 
. 15000 

.00540 

.00540 

160.0 

70.0 


.14932 
. 14999 

.00540 

.00540 

160.4 

70.3 


.451 % 
.007 % 

.000001 Hz 
.000001 Hz 

.45 

.31 


.069 

.096 


.79 

.78 


.500 

59.0 


Pendulous data, amp, freq, P h “® , 

Equivalent U t, V deflections (meters) . 

lSe ?-: e dat;- 2000 ? sigma value (deg/sec, : -000280 

g t op° indexes tart index, 8 number of point - 593 


.03125 

58.7 


50.0 



REV LABEL 
0 
1 
2 

3 

4 


ME FOR THIS 

DATA WINDOW 

IS 

• 

• 

. 303104E+03 

[ LAST POINT 

IN BUFFER IS 

# 

• 


. 102298E+04 

POLARITY 

RATE (D/S) 

# 

OF ROTATIONS START TIME 

POSITIVE 

1.9437 


5.9 

1190.51 

POSITIVE 

1.9437 


5.9 

1375.72 

POSITIVE 

1.9437 


5.9 

1560.93 

POSITIVE 

1.9437 


5.9 

1746 . 14 

POSITIVE 

1.9437 


5.9 

1931.35 


Amp x 
Amp y 

Freq x 
Freq y 

Phase x 
Phase y 


INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

. 15000 
. 15000 

. 14946 
. 14998 

.357 % 
.016 % 


.00540 

.00540 

.00540 

.00540 

.000001 Hz 
.000002 Hz 

i 

160.0 

70.0 

160.3 

70.5 

.27 

.46 

.055 

.150 

data, amp, 

freq, phase 

— * 
• 

.500 

pa a c a 


PHASE 0 T=15 MIN 
degs 


.53 

1.19 


50.0 


lseed = -130000 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point = 593 


.000280 
1 593 



midpoint time for this data 
time tag on last point in buffer is • 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9436 
1.9436 
1.9436 
1.9436 
1.9436 


# 


: . 303 104E+03 

. 102298E+04 

OF ROTATIONS START TIME 
5 9 1190.60 

5 ' 9 1375.82 

5 ’ 9 1561.05 

5 * 9 1746.27 

5*9 1931.50 



INPUT 

OUTPUT 

ERROR 

Amp x 
Amp y 

. 15000 
. 15000 

. 14980 
. 15020 

.136 % 
-.134 % 

Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000003 Hz 
.000000 Hz 

Phase x 
Phase y 

160.0 

70.0 

160.5 

70.1 

.46 
. 15 


PHASE GRAD 
degs/cycle 


. 181 
.025 


PHASE 6 T=15 MIN 
degs 


1.33 

.27 


Pendulous data, amp, freq, phase _• 

Equivalent U & V deflections (meters) -. 

.Led = -140000 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point 5 


.500 

59.0 


.03125 
58 . 9 


50.0 


.000280 
1 593 



MIDPOINT 

TIME FOR THIS 

DATA WINDOW 

TIME TAG 

ON LAST POINT 

IN BUFFER IS 

REV LABEL POLARITY 

RATE (D/ S ) 

0 

POSITIVE 

1.9439 

1 

POSITIVE 

1.9439 

2 

POSITIVE 

1.9439 

3 

POSITIVE 

1.9439 

4 

POSITIVE 

1.9439 


: . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 

5.9 1190.45 

5.9 1375.65 

5.9 1560.85 

5.9 1746.04 

5.9 1931.24 



INPUT 

OUTPUT 

ERROR 


PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

. 15000 
. 15000 

. 14995 
. 15008 

.033 % 
-.053 % 



Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000000 

.000001 

Hz 

Hz 

t 

Phase x 
Phase y 

160.0 

70.0 

160.4 

70.3 

.37 

.28 


.019 

.075 


PHASE @ T=15 Mil 
degs 


.46 

.64 


Pendulous data, amp, freq, phase : 

Equivalent U & V deflections (meters) =: 
iseed = -150000 . . 

Noise data: 1 sigma value (deg/sec) 

Stop index, start index, & number of point 


.500 

59.0 


.03125 

58.9 


50.0 


= 593 


000280 
1 593 


MIDPOINT 

TIME FOR THIS 

DATA WINDOW 

TIME TAG 

ON LAST POINT 

IN BUFFER IS 

REV LABEL POLARITY 

RATE (D/S) 

0 

POSITIVE 

1.9439 

1 

POSITIVE 

1.9439 

2 

POSITIVE 

1.9439 

3 

POSITIVE 

1.9439 

4 

POSITIVE 

1.9439 


s . 303104E+03 

. 102298E+04 

OF ROTATIONS START TIME 


5.9 

5.9 

5.9 

5.9 

5.9 


1190.44 

1375.64 

1560.83 

1746.03 

1931.23 



INPUT 

OUTPUT 

ERROR 


PHASE GRAD 
degs/cycle 

PHASE @ T=15 MIN 
degs 

Amp x 
Amp y 

. 15000 
. 15000 

. 14975 
. 14999 

.165 % 
.009 % 




Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000000 

.000001 

Hz 

Hz 



~ Phase x 
Phase y 

160.0 

70.0 

160.3 

70.4 

.27 

.40 


.022 

.071 

.37 

.74 


— Pendulous data# amp, freq, phase 

Equivalent U & V deflections (meters) =: 58. 

iseed = -160000 

_ Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point = 593 


.03125 

58.9 

000280 
1 593 


50.0 



- MIDPOINT TIME FOR THIS DATA WINDOW IS 
TIME TAG ON LAST POINT IN BUFFER IS 

A rt « mT-t / Pk / C \ 


. 303104E+03 
. 102298E+04 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 

POSITIVE 


RATE (D/S) 
1.9439 
1.9439 
1.9439 
1.9439 
1.9439 


# OF ROTATIONS 
5.9 
5.9 
5.9 
5.9 
5.9 


START TIME 
1190.44 
1375.64 
1560.84 
1746.04 
1931.24 



INPUT 

OUTPUT 

ERROR 

PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

. 15000 
. 15000 

. 14961 
. 15005 

.260 % 
-.032 % 


Freq x 
Freq Y 

.00540 

.00540 

.00540 

.00540 

.000000 Hz 
.000001 Hz 


Phase x 
Phase y 

160.0 

70.0 

160.3 

70.2 

.33 

.20 

.003 

.056 


PHASE 0 T=15 MIN 
degs 


.34 

.47 


Pendulous data, amp, freq, phase • 

Equivalent U & V deflections (meters) 
iseed = -170000 

Noise data: 1 sigma value (deg/sec) 

Stop index, start index, & number of point 


.500 

59.0 


.03125 

58.8 


50.0 


= 593 


000280 
1 593 



MTDPOINT time for this data window is 

TIME^TAGON^LAST^POINT » ™ « ' 

positive 

POSITIVE 
POSITIVE 
POSITIVE 
POSITIVE 


REV LABEL 
0 
1 
2 

3 

4 


RATE(D/S) 
1.9436 
1.9436 
1.9436 
1.9436 
1.9436 


.303104E+03 
. 102298E+04 
# OF ROTATIONS START TIME 

1190.60 
1375.82 
1561.04 
1746.27 
1931.49 


5.9 

5.9 

5.9 

5.9 

5.9 


INPUT 

OUTPUT 

ERROR 

. 15000 
.15000 

.14971 
. 14962 

.195 % 
.255 % 

.00540 

.00540 

.00540 

.00540 

.000002 
. 000000 

160.0 

70.0 

160.4 
70 . 3 

.44 

.33 


Amp x 
Amp y 

Freq x 
Freq Y 

phase x 
Phase y 

pendulous data, amp, freq, . 

Equivalent U & V deflections (meters) . 

iseed ” j. t I^ 0OO °i sigma value (deg/sec) : 
Stop 0 indexes tart InUl 4 numter of point 


PHASE GRAD 
degs/cycle 


PHASE § T=15 Mil 
degs 


.148 

.006 


1.16 

.35 


.500 

58.8 


.03125 

58.8 


50.0 


= 593 


.000280 
1 593 


REV LABEL 
0 
1 
2 

3 

4 


ME FOR THIS 

DATA WINDOW 

IS 

• 

• 

. 303104E+03 

[ LAST POINT 

IN BUFFER IS 

• 

• 


. 102298E+04 

POLARITY 

RATE (D/S) 

# 

OF ROTATIONS START TIME 

POSITIVE 

1.9438 


5.9 

1190.46 

POSITIVE 

1.9438 


5.9 

1375.67 

POSITIVE 

1.9438 


5.9 

1560.87 

POSITIVE 

1.9438 


5.9 

1746.07 

POSITIVE 

1.9438 


5.9 

1931.27 



INPUT 

OUTPUT 

ERROR 


PHASE GRAD 
degs/cycle 

Amp x 
Amp y 

. 15000 
. 15000 

. 14941 
. 15026 

.391 % 
-.172 % 



Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000001 

.000000 

Hz 

Hz 


Phase x 
Phase y 

160.0 

70.0 

160.6 

70.2 

.55 

.23 


.087 

.020 


PHASE 6 T=15 MIN 
degs 


.98 

.33 


Pendulous data, amp, freq, phase 
Equivalent U & V deflections (meters) 
iseed = -190000 

Noise data: 1 sigma value (deg/sec) : 

Stop index, start index, & number of point — 


.500 .03125 

59.1 58.7 

.000280 
593 1 593 


50.0 


midpoint time for this data window is 
time°tag on last point in buffer is . 


REV LABEL 
0 
1 
2 

3 

4 


POLARITY 

positive 

positive 

positive 

positive 

positive 


RATE (D/S) 
1.9431 
1.9431 
1.9431 
1.9431 
1.9431 


# OF ROTATIONS 
5.9 
5.9 
5.9 
5.9 
5.9 


. 303104E+03 
. 102298E+04 


START TIME 
1190.84 
1376 .11 
1561.38 
1746.64 
1931.91 



INPUT 

OUTPUT 

ERROR 

Amp x 
Amp y 

. 15000 
. 15000 

. 14966 
. 14979 

.229 % 
.139 % 

Freq x 
Freq y 

.00540 

.00540 

.00540 

.00540 

.000002 Hz 
.000003 Hz 

Phase x 
Phase y 

160.0 

70.0 

160.5 

70.6 

.54 

.58 


Pendulous data/ amp/ freq, phase 
Equivalent U & V deflections (meters) 

1SS Noise dataf 000 ? sigma value (deg/sec) : 
Stop index, start index, & number of point 


PHASE GRAD 
degs/cycle 


PHASE @ T=15 MIN 
degs 


. 149 
.169 


1.26 

1.40 


.500 

58.9 


.03125 

58.8 


50.0 


593 


.000280 
1 593 



INPUT OUTPUT 


ERROR 


PHASE GRAD 
degs/cycle 


PHASE § T=15 MIN 
degs 


Amp x 
Amp Y 


. 15000 
. 15000 


Freq x 
Freq y 


.00540 

.00540 


Phase x 
Phase y 


160.0 

70.0 


. 14966 
. 15003 


.226 % 
-.017 % 


.00540 

.00540 


.000001 Hz 
.000001 Hz 


160.4 

70.3 


.43 

.26 


freq, phase 


Pendulous data, /^^tiinrimeters) =: 

S top° indexes tart Intel * number of poxnt - 


.074 

.069 

, 500 .03125 

59.0 58.8 

.000280 
593 1 593- 


.79 

.60 


50.0 


INPUT OUTPUT 


ERROR 


PHASE GRAD 
degs/cycle 


PHASE @ T=15 MIN 
degs 


Amp x 
Amp y 


. 15000 
. 15000 


Freq x 
Freq y 


.00540 

.00540 


Phase x 
Phase y 


160.0 

70.0 


.14965 
. 15001 


.232 % 
-.006 % 


.00540 

.00540 


.000001 Hz 
.000001 Hz 


160.4 

70.3 


.41 

.30 


pendulous data , ^ i 

EquiYalent U & d “ value (deg/sec) s 

Stop 0 indexes tart Indel « number of point - 


.075 

.072 

.500 .03125 

59.0 58.8 

.000280 
593 1 593 


.77 

.64 


50.0 



APPENDIX 2 . C 


Programs 

( 1 ) 

( 2 ) 

( 3 ) 


for ECR Testing 
SIMREAD - Read Simulation 
SIMTEST - Compare Observer 
BTBUSO - Compare Observer 


Data Files 

Output to Simulation Input 
Output to Model Signal Input 



n o ri n onooonortooo 


The program SIMREAD . FOR reads a simulation input file and 5 re *^® s J" 1 ?® 
output file IDFTXY.DAT for the UNOMSC.FOR program (frequency domain skip 

rope observer) • 


1 


2 


real x,y, tlength, time 
integer amcsmode 

open (10, file = ' TRUTH . DAT ' , status = 'old ) 
open (11, file = 'IDFTXY.DAT' , status = 'unknown ) 

read^lO, * , end =2) time, x, y, tlength, amcsmode, dl, d2 
write(ll,*) time, x, y, tlength, amcsmode 
i = i + 1 


go to 1 
continue 
If = i - 1 

print*, 'read ',lf,' lines of data' 
close ( 10 , status = 'keep') 
close (11, status = 'keep') 
end 


The program SIMTEST. FOR compares the output of UNOMSC.FOR to the ongi 
nal simulation input data file. 


PROGRAM TO TEST FILTER AGAINST SIMULATION DATA 

PHASE ANGLE AT TIME T IS DEFINED AS ARC TAN (V TERM/ U TERM) 
PHASE ERROR IS DEFINED AS 'TRUTH' - 'MODEL'. 

MAGNITUDE CALCULATED AS SQRT (U**2 + V**2 ) 

MAGNITUDE ERROR EXPRESSED AS PERCENT TERM BY /mr , rTm „ % 

ERROR. = 100% * ( MAG (TRUTH) - MAG (MODEL) )/ MAG (TRUTH ) 

FINAL ERRORS ARE EXPRESSED AS RMS OVER ALL TIME POINTS 


NAME IS 'SIMTEST' 

INTEGER NDIM 
REAL* 4 DFR 

PARAMETER (NDIM=4000, DFR=57.2958 ) 


REAL* 4 UM(NDIM) , VM (NDIM) , UT(NDIM) , VT (NDIM) 
REAL* 4 RE (NDIM), RM (NDIM) , RT (NDIM) 

REAL* 4 PE (NDIM), PM (NDIM), PT (NDIM) 


OPEN INPUT FILES - 'F TXYUV.DAT' IS OUTPUT FROM 'UNOMSC.FOR' 

?MUST SET ODF TIME TO _ .TRUE. TO WRITE THIS FILE ! ! ) AND ' TRUTH . DAT ' 
IS THE SIMULATION DATA FILE (COPY THE APPROPRIATE SIMULATION DATA 
FILE INTO ' TRUTH . DAT ' BEFORE RUNNING ' SIMTEST. FOR' ) 


OPEN (11, FILE = 
OPEN (12, FILE = 


'F TXYUV.DAT', STATUS = 'OLD') 
' TRUTH . DAT ' , STATUS = 'OLD') 


1 



CALCULATE THE TIME SAMPLING RATES AND 


READ THE INPUT FILES, AND 
LENGTH OF EACH 


READ (11, *, 
READ (11, *, 


END=3) TIME1, 
END=3) TIME2 , 


Dl, D2 , UM ( 1) , 
Dl, D2 , UM(2), 


VM ( 1) 
VM (2) 


READ (11, *, END=3) TIME, Dl, D2 , UM(I), VM(I) 
1 = 1 + 1 
GO TO 2 
LM = I “ 1 

DT M = TIME2 - TIME1 


READ 

(12, 

*, 

END=7 ) 

READ 
1 = 3 

(12, 

*, 

END=7 ) 

READ 

(12, 

*, 

END=7) 


1 = 1 + 1 
GO TO 6 
LT = I “ 1 
dt T = T1 - TO 


TO, Dl, D2 , D3 
Tl, Dl, D2 , D3 

T, Dl, D2 , D3 , 


D4 , UT(1) , VT(1) 
D4 , UT (2 ) , VT (2 ) 

D4 , UT (I) , VT (I) 

} 


LC = MINO ( LM, LT) 
RLC = FLOAT (LC) 


PRINT * , ' 
PRINT * , ' 
PRINT * , ' 
PRINT * , ' 
PRINT * , ' 


yiODEL DATA HAS ' , LM, ' TIME 
TRUTH DATA HAS ', LT, ' TIME POINTS 
WILL USE ' , LC, ' TIME POINTS 
DT FOR MODEL DATA IS : ' , DT_M 
DT FOR TRUTH DATA IS : ' , DT_T 


9 

9 


WRITE (7,*) 
WRITE ( 7 , * ) 
WRITE (7 , *) 
WRITE (7,*) 
WRITE (7 , *) 


' MODEL DATA HAS ' , 

' MODEL DATA HAS ' , 

' WILL USE ' , LC, ' 
' DT FOR MODEL DATA 
/ DT FOR TRUTH DATA 


LM, ' TIME POINTS 
LT, ' TIME POINTS 
TIME POINTS ' 

IS : ' , DT_M 
IS : ', DT_T 


9 

9 


CALCULATE THE 
DO K=l, LC 

RM(K)= SQRT( 
RT (K) = 

PM (K) = 
PT(K)= 

RE (K ) — 

END DO 


OVERALL RMS ERRORS IN THE MAGNITUDE AND PHASE 


SQRT( 
DFR * 
DFR * 
100 . 


UM (K) **2 + VM (K) **2 ) 
UT ( K ) * * 2 + VT (K) **2 ) 
AT AN 2 ( VM(K) , UM(K) ) 
AT AN 2 ( VT (K) , UT(K) ) 
* (1. - RM(K)/RT(K) ) 


CALL SATAN (PM, LC, 1) 
CALL SATAN (PT, LC, 1) 


SSQP = 0. 

SSQR = 0. 

DO K=1 , LC 

2 



PE (K) = PM(K) - PT(K) 
SSQP = SSQP + PE (K) **2 
SSQR = SSQR + RE (K) **2 
END DO 


PHASE = SQRT ( SSQP/RLC ) 

RMAG = SQRT ( SSQR/RLC ) 

PRINT * , ' OVERALL RMS MAGNITUDE ERROR 
PRINT * , ' OVERALL RMS PHASE ERROR ~ 


RMAG, ' %' 

PHASE, ' DEGREES' 


WRITE (7,*) ' OVERALL RMS MAGNITUDE ERROR 

WRITE (7,*) ' OVERALL RMS PHASE ERROR 


RMAG, ' %' 
PHASE,' DEGREES' 


C 


CLOSE FILES 

CLOSE (11, STATUS = 'KEEP') 
CLOSE (12, STATUS - 'KEEP') 


END 



C 


C 

C 

5 


C 


SUBROUTINE SATAN (A, N, TYPE) 


’MART ATAN PROGRAM - REMOVES THE 2 PI JUMPS AT -+ 180. 

THE PROGRAM IS A POST PROCESSOR. IT COULD BE 

iLTEREyO ; N degrees , both INPUT AND OUTPUT. 

[F TYPE NE 1, THEN RADIANS WILL BE USED. 


REAL* 4 A ( 1) , CHECK, ADD, PI 
INTEGER TYPE, N , JK, JKO 
LOGICAL FLAG 

PARAMETER (PI = 3.1415926) 


IF (TYPE .EQ. 1) THEN 
CHECK = 355.0 
ADD = 360.0 
ELSE 

CHECK =6.19 
ADD = 2.0 * PI 
END IF 


JKO = 2 

FLAG = . TRUE . 
JK = JKO 


, GT . CHECK) THEN 


DO K = JK, N 

IF (ABS (A (K) - A (K-l) ) 

A(K) = A(K) + SIGN (ADD, A(K-l)) 

IF (FLAG) THEN 
FLAG = .FALSE. 

JKO = K 


3 



END IF 
END IF 
END DO 


IF (.NOT. FLAG) GO TO 5 

RETURN 

END 

The program VTBUSO.FOR compares the results 
signals generated by CREATE. FOR. 


of UNOMSC . FOR .to the model 


c 


c 

c 

c 


c 


open(3 , f ile='F REC.DAT' , status='old ) 
open ( 4 , f i le= ' F_RECORD . DAT ' , sta tus= old ' ) 
open ( 7 , f i le= ' T_REPORT . DAT ' , s tatus= unknown ) 


read file from signal generating program (CREATE. FOR) 

read (3,*) knoise, tensigma, iseed 

read (3 , *) axin, fxin, pxin 

read (3 , *) ayin, fyin, pyin 

read (3 , *) apend, fpend, ppend 


output 


make sure angles are bounded by pi. 
subroutine a j sign puts angles in 


-pi to + pi range. 


pxin = a j sign (pxin) 
pyin = ajsign (pyin) 

frqin = 0.5 * (fxin + fyin) 
period = 1.0/frqin 
onesigma = tensigma/ 10. 0 

read file from filter output 

read (4,*) tm 

read (4 , *) lb, le, leb 

read ( 4 , * ) to, tf, dt, tl, wk, psign 

read (4 , *) ax, fx, px, ay, fy, py 

read (4 , *) fa, ta, u, v, fl, fh 

read (4, *) avgx, avgy 


px = 57 . 296*px 
py = 57 . 296*py 


epx = abs ( abs( pxin ) - abs ( px) ) 
epy = abs ( abs( pyin ) - abs( py) ) 
if (axin .ne. 0.0) then 

eax = 100. *(1. - ax/axin) 
else 

eax = 357 . 
end if 


4 



if (ay in .ne. 0.0) 
eay — 100 . * ( 1 . - 
else 

eay = 357. 
endif 


then 

ay/ayin) 


c 

c 


9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
20 
22 


(900. /period) 


ax, 

ay, 


eax 

eay 


pgx = 360 . *abs (frqin - fx) *period 
pay = 360 . *abs (frqin - fy) *period 
phase_errx = epx + pgx * (900 . /period) 
phase_erry = epy + pgy * 
close ( 3 ,status='keep' ) 
close ( 4 , status= ' keep ' ) 
write (7 ,9) 
write (7, 9) 
write (7, 9) 
write(7 ,9) 
write(7 , 10) 
write (7 , 11) axin, 
write (7 ,12) ay in, 
write (7 ,9) 
write (7 , 13) 
write (7, 14) 
write (7, 9) 
write (7 , 15) pxin, 
write (7, 16) pyin, 
write (7, 9) 
write (7 , 17) 
write (7 , 18) 
if (knoise 
write (7 , *) 
write(7 , 20) 
else 

write (7 , *) 
endif 

write (7 , 23) le, lb, leb 


frqin, fx, abs (frqin-fx) 
frqin, fy, abs (frqin-fy) 


px, 

py, 


epx, 

epy, 


pgx, 

pgy, 


phase_errx 

phase_erry 


apend, fpend, ppend 
u, v 

eq. 1) then 

iseed = ' , iseed 
onesigma 

No noise for this run.' 


format (12X, ' INPUT ' , 6X , ' OUTPUT ' , 8X , ' ERROR' , 6X, ' PHASE GRAD ' 

$ 6X 'PHASE @ T=15 MIN' , / 48 x, 'degs/cycle' , 12x, degs ,/) 
format (2X, 'Amp x' , f 10 . 5 , f 12 . 5 , f 12 . 3 , ' %') 

format (2X, 'Amp y ' , f 10 . 5 , f 12 . 5 , f 12 . 3 , ' % ) 

format ( 2 x, 'Freq x' , f 9. 5, f 12 . 5, f 14. 6, Hz ) 
format ( 2 x, 'Freq y ' , f 9 . 5 , f 12 . 5, f 14 . 6, Hz ) 
format (2x, 'Phase x' , f 8 . 1 , f 12 . 1, f 12 . 2 , 6x, f 8 . 3 , 10x, f 8 . 2) 
format(2x, 'Phase y' , f 8 . 1 , f 12 . 1 , f 12 . 2 , 6x, f 8 . 3 , lOx, f . ) 
format(2x, 'Pendulous data, amp, freq, phase -. 

^format*( 2 x^ Equivalent U 6 V deflections (»®ters)=: ; , 2f 9 . 1) 
fonnat(5x, 'Noise data: 1 sigma value (deg/sec) • <f 14 - 6 > 

format ( 2 x, 'Means in X and Y axes were :',2fl2.7, and 
$ , ' were removed') 
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, & number of point -',314) 


23 format ( 2 x, 'Stop index, start index 

write(7 ,9) 
write (7 ,9) 

close (7, status= ' keep ' ) 

call exit 

end 

real function ajsign (x) 
if ( abs(x) .gt. 180. ) then 
ajsign — x - sign(360., x) 
else 

ajsign = x 
end if 
return 
end 
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APPENDIX 2 . D 


Programs 

( 1 ) 

( 2 ) 

(3) 

(4) 


for Systematic Testing 

NO NOISE - Noise-free Test Cases, Station 2 
NoTse - Noisy Test Cases, Station 2 - page 
LIBRATION - Noise-free Tests at Station 1 
LIB NOISE - Noisy Tests at Station 1 - page 


- page 1 

page 16 
23 



_ mo for uses the model signal creation (without noise) algorithm 

Britan and the WORK, HANN, MEAN, F0UR1, and LSCF sub- 
c routines from the UNOMSC.FOR program. NO_NOISE.FOR systematic y 
_C through the phase pairings for a given frequency. 


REAL X (3000 ) 

COMMON / PARI / LB , LE , NFT , NPNT , NDEG 

COMMON/ PAR2 / DT , PI , DF , PID 

dt = 1.024 

sd = 2 . 8e-03 

read*, fix, iseed 

lb=l 


npnt = 7 
ndeg = 3 

le = int(l/ (f lx*dt) +0 . 5) 
if (le.eq. 2* (le/2) ) le = le+1 
le = le*3 + lb - 1 


NFT=8192 
LEB=LE-LB+1 
PI=4 . 0 * ATAN (1.0) 

DF=1 • 0/ (NFT*DT) 

PID=180. 0/PI 
do m = 1,37 

phlx = (m-19) *10 
do n = 1,37 

ph2x = (n-19) *10 
CALL CREATE (X , A1X , FIX , PH1X, ph2x) 
CALL WORK (X, ampx, f rx, phx, aerror , 
write ( 11 , * ) phlx , ph2x , aerror 
write (12 , *) phlx,ph2x, f error 
write ( 13 , * ) phlx , ph2x , perror 
end do 
end do 


f err or , perror , A1X , FIX , PH1X) 


END 


SUBROUTINE CREATE ( xt , A1X , FIX , PH1X , ph2x) 

REAL xt (3000) 

dt = 1.024 

iper=1000 

alx=0. 02 

a0x=0 . 065 

a2x=0 . 5 

PI = 4 . 0*ATAN (1.0) 
convrt = pi/180.0 
phi lx = phlx*convrt 
phi2x = ph 2 x*convrt 
TPIDT = 2 . 0*PI*DT 
F2X = 0.03125 
xtheta = tpidt*f lx 
xphi = tpidt*f2x 
DO I = 1 , IPER 
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II 


xt(I) = AOX + A1X*C0S (xtheta*Il+PHilX) +A2X*COS (xphi*Il+PHi2X) 

END DO 
RETURN 
END 


C 

C 

C 

C 

C 

C 

C 


SUBROUTINE WORK CALCULATES THE AMPLITUDE, PHASE, AND FREQUENCY 

OF THE DATA THE FOURIER TRANSFORM SUBROUTINE FOUR1 IS 

oa T TFD BY SUBROUTINE WORK. WORK RETURNS TO THE MAIN PROGRAM 

Values of the amplitude, phase, and frequency as well as 

THE TIME INDEX WHERE THE MAXIMUM VALUE OCCURS. 

THIS IS BASED ON MODEL OF COS (WT+PHASE) . 

SUBROUTINE WORK (ANG , amp , f req , phase , Aerror , Perror , Ferror , A1 , FI , PHI1) 

PARAMETER D ™DIM-3000 , NCDIM=8200, NFT-8192, NPNT-7) 

DIMENSION AUX(NDIM) ,ANG(1) 


REAL* 4 XFREQ(7), PHIMAG (7 ) , PHREAL (7 ) 
COMPLEX AWO(NCDIM) 


C 

c 

c 

c 

c 


COMMON / PARI / LB , LE 
COMMON/ PAR2 / DT , PI , DF , PID 


NTB1=LE-LB+1 

NTB1 IS FORCED TO BE ODD 

HANN WINDOW ROUTINE USES 


IN MAIN PROGRAM. 

ODD NUMBER OF POINTS TO TAPER. 


LOAD INPUT DATA FROM ANG (I) INTO ARRAY AUX(J). 
LB 1=1 “LB 
DO I=LB , LE 
IL=I+LB1 
AUX ( IL) =ANG ( I ) 

END DO 


APPLY WINDOW FUNCTION TO TIME SEQUENCE 
CALL HANN (NTB1 , AUX , BIAS) 


C 

C 

C 

C 


AKE COMPLEX NUMBER AWO(I) FROM R ^fiF ^ 5 

ZERO IMAGINARY VALUE (AUX (I) IS THE REAL VALUE). 


BY USING 


DO 1=1 , NTB1 

AWO ( I ) =CMPLX ( AUX ( I ) ,0.) 
END DO 


NOW PAD THE DATA STREAM WITH ZEROS OUT TO AWO (8192) . 
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DO I=NTB1+1,NFT 

AWO ( I ) = CMPLX ( 0 . , 0 . ) 
END DO 


C SUBROUTINE FOUR1 DOES THE FOURIER TRANSFORM USING A FFT METHOD 

C CALL FOUR1 (AWO , NFT , 1 ) 

c LOOP TO FIND THE MAXIMUM MODULUS VALUE OF THE FOURIER TRANSFORM 
ampMAX=0.0 

istart = int ( (f 1-0. 001) /df)+l 
lend = int ( (f 1+0 . 001 ) /df ) +1 
do i = istart, iend 
f5T = (i — 1) 

if (cabs(awo(i) ) . gt . ampmax ) then 
ampmax - cabs(awo(i)) 
kf = i 
freq = fr 
end if 
end do 



C 


Sm™l Is*™ degree°and "S?nt B s Y wi“ S Ie S Ss^ E in P 2™ve M fit: 

3 ^MAGNITUDE OF TRANSFORM (SQRT (REAL**2 + IMAG**2)) 

REAL PART 
IMAGINARY PART 

CENTER OF DATA SET IS THE FREQUENCY POINT WHERE MAX WAS FOUND. 

DO I = 1 , NPNT 

J = KF- ( (NPNT+1) /2 . 0) +1 
XFREQ(I) = CABS ( AWO ( J) ) 

PHIMAG(I) = AIMAG ( AWO ( J) ) 

PHREAL ( I ) = REAL ( AWO ( J ) ) 

END DO 


>0 CURVE FIT ON THE MODULUS OF THE FOURIER TRANSFORM 
•ATL TO LSCF WITH OPTION 1 DOES 2 THINGS. »-,/-> -rvim 

CURVE FITS AND COMPUTES TRUE MAX FREQUENCY POINT. 


CALL LSCF (FQ_P0, ampMAX, XFREQ , 1) 


;ALL LSCF (FQ_P0, PHASEI, PHIMAG, 
2ALL LSCF (FQ_P0, PHASER, PHREAL, 
?REQ = FREQ + DF * FQ_P0 


2 ) 

2 ) 


SCALING OF TRANSFORMED DATA IS PERFORMED TO GIVE OUTPUTS IN 
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c 


- c 
c 
c 

_ c 
c 
c 
c 

“ c 
c 
c 

- c 
c 
c 


DEGS/SEC AND REPRESENT ACTUAL RATE DATA. 
SCALE = 4 . O/FLOAT (NTB1-1) 

AMP = SCALE * ampMAX 

PHASE = -ATAN2 (PHASEI , PHASER) *pid 


THE FORWARD FOURIER TRANSFORM IS USUALLY ^FINED 

i , and**we°use 

"--IS sr 

S& c SST^JS5JS5J 5 S«^SS'J5^^ .. 


FERROR = (ABS(FREQ-Fl) / 0.005) * 100 

kERROR = (ABS(AMP-A1)/A1) * 100 

if^abs^philfleq^lSOJ^perror^^abstabs (phase) -abs (phil) ) / 360. 0) *100. 


return 

END 


0 


SUBROUTINE HANN (LA, All, BIAS) 


TAPER IS RAISED COSINE CURVE. 

MEAN IS COMPUTED AND REMOVED FROM INPUT SIGNAL 
PARAMETER (PI=3 . 1415926) 

REAL All (LA) , HW(3000) 

ITM = (LA-1) / 2 
RM = FLOAT (ITM) 

DO IT= -ITM, ITM 

I = 1 + IT + ITM % „ 

HW(I)=0.5*(1.0 + COS (PI*FLOAT (IT) /RM ) ) 

All (I) =A11 (I) *HW(I) 


END DO 

toSe^EaS 1 ^ TWICE P COMPUTED N VALUE BECAUSE HANN WINDOW 
REDUCES^VALUE BY FACTOR OF 2 . (I.E. MEAN OF WINDOW IS 


0.5) 


CALL MEAN (LA, All, BIAS) 

BIAS = BIAS * 2.0 * LA/ (LA-1) 

DO I = 1 , LA 

All ( I) = All (I) - BIAS * HW(I) 
END DO 
RETURN 
END 
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SUBROUTINE MEAN ( LA , A2 2 , SA) 

THIS ROUTINE COMPUTES THE DC TERM OF THE DATA STREAM. 
MEAN IS NOT REMOVED, BUT ONLY COMPUTED. 

REAL A22 (LA) 

SA = 0. 

DO 1=1, LA 

SA=SA+A22 (I) 

END DO 

SA=SA/ FLOAT (LA) 

RETURN 

END 


SUBROUTINE FOUR1 (DATA,NN, ISIGN) 

THIS ROUTINE DOES THE FOURIER TRANSFORM USING A FFT METHOD. 

REAL* 8 WR, WI ,WPR,WPI ,WTEMP, THETA 
DIMENSION DATA ( * ) 

N=2*NN 
J=1 

DO 11 1=1, N, 2 
IF ( J. GT . I) THEN 
TEMPR=DATA ( J ) 

TEMPI=DATA ( J+l ) 

DATA(J)=DATA(I) 

DATA (J+l) =DATA ( 1+1) 

DATA(I)=TEMPR 
DATA ( I + 1 ) =TEMP I 
ENDIF 
M=N/2 

IF ( (M.GE.2) .AND. (J.GT.M) ) THEN 
J=J-M 
M=M/2 
GO TO 1 
ENDIF 
J=J+M 
11 CONTINUE 

MMAX=2 

2 IF (N.GT.MMAX) THEN 

ISTEP=2*MMAX 

THETA=6 . 283 1853 07 17 95 9D0/ (ISIGN*MMAX) 
WPR=-2.D0*DSIN(0.5D0*THETA) **2 
WPI=DSIN (THETA) 

WR=1.D0 
WI=0 . DO 

DO 13 M=1 ,MMAX , 2 
DO 12 I=M, N , I STEP 

J=I+MMAX , , , . 

TEMPR=SNGL (WR) *DATA ( J) -SNGL(WI) *DATA(J+1) 
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12 


13 


TEMPI=SNGL(WR) *DATA( J+l) +SNGL(WI) *DATA( J) 

DATA(J)=DATA(I) -tempr 

DATA( J+l) =DATA(I+1) -TEMPI 

DATA ( I ) =DATA ( I ) +TEMPR 

DATA ( 1+1) =DATA (1+1) +TEMPI 


CONTINUE 

WTEMP=WR 

WR=WR*WPR-WI *WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 


CONTINUE 

mmax=istep 
GO TO 2 
ENDIF 
RETURN 
END 



f5^A S MD°MGrL D ?0LYN0MiL SQ ™ E Da5Hs ASSraED 7 TO°BflAHPLED 
[? R integ^l in?Ir»ls m.y scaling must be done outside 
THIS SUBROUTINE. THE 7 POINTS ARE . 

the”polynomial 1 is 0 f(P) -'a + b *p + J*p*p. 

THE MAX OCCURS AT P - PO - "B/ (2 C) . DATA str£AM) 

this^delt A C is^referenced to midpoint frequency OF 7 POINTS. 
™e MAX VALUE IS F(PO) = A - (B*B) / (4*C) . 


SUBROUTINE LSCF (PO, FMAX, U_IN, IOPT) 



ON ENTRY I 

U IN IS INPUT ORDINATE VALUES* (7 ) 

l5P T " S p?Nn T ?0 WERE mX 2 OCCTOS S RLUS COMPUTE MAX VALUE. 

2 • COMPUTE VALUE OF POLYNOMIAL AT SPECIFIED FREQUNCYPO. 

IF IOPT=°V THEN PO IS FREQUENCY POINT TO EVALUATE POLYNOMIAL. 


PO IS VALUE OF P WHERE MAX PEAK OCCURS; 

THIS IS WRT CENTER POINT OF DATA. 

FMAX IS VALUE OF FUNCTION AT P=P °* .. 

*********************************** 

REAL* 4 U_IN ( * ) 

r******* ^** ****** ******************* ******* 

wtrqt step is to do least squares* 

all COEFFICIENTS HAVE BEEN P * E_C f^!£ED by 
'A' IS -8, 12, 24, 28, 24, 12, -8 DIVIDED BY 
,t>/ ts —9 -6. -3, 0, 3, 6, 9 DIVIDED BY 

IS 5 ,'o. -3. -«. -3. 0, 5 DIVIDED BY 


84 

84 

84 


US17 = U IN(1) + U_IN(7) 
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c 


c 


c 


US35 = U_IN ( 3 ) + U_IN ( 5 ) 

A = C0F1 

COFl = -8 . *US17 + 12. * (U_IN(2)+U_IN(6) ) + 
l 24 . *US3 5 + 28 . *U_IN (4) 

B = COF2 

COF2 = 9.*(-U IN(1) +U_IN(7) ) + 
f 6.*(-U_IN(2)+U_IN(6) ) + 

S 3 . * ( -U_IN ( 3 ) +U_IN ( 5 ) ) 


c== qqp2 

COF3 = 5 . *US17 -3 . *US3 5 - 4.*U_IN(4) 

IF (ABS(C0F3) .LT.1.0E-08) THEN 

PRINT *, ********** WARNING ********* 
PRINT* , ' CONSTANT VALUE EQUALS ZERO 
PRINT* , ' FREQUENCY VALUE . ' 


- CANNOT COMPUTE 


return 

ENDIF 


A MAXIMUM' 


DID NOT DIVIDE BY 84 YET. DO SO FOR MAX PART BUT NOT PO 


COMPUTE PO, VALUE OF F( 
COMPUTE FUNCTION AT PO; 
IF (IOPT .EQ. 1) THEN 
PO=-0.5*COF2/COF3 
FMAX = (COFl + 0.5 * 


P) WHERE A+B*P+C*P*P = 0 
A+B*P0+C*P0*P0 = A-B**2/4C 


PO * COF2 ) /84 . 


FMAX = (COFl + P0*( COF2 + PO*COF3 ) )/84. 

END IF 

RETURN 

END 


C 

c 

c 

c 


ror?L“Sm'i'E 0 FOR S p?og«m S and a thf , ^FOTOl^'and^lcF^sub- 

outines from the UNOMSC.FOR program. NOISE. FOR systematically runs roug 
he phase pairings for a given frequency. 


REAL X (3000) ,xinit(3000) ,ax(100) ,fx(100) ,px(100) ,axe(100) ,fxe(100) , 
# pxe(100) 

EXTERNAL rani 

COMMON / PARI / LB , LE , NFT , NPNT , NDEG 

COMMON /PAR2 /DT , PI , DF , PID 

dt = 1.024 

sd = 2 . 8e-03 

read*, fix, iseed 

lb=l 

npnt = 7 
ndeg = 3 

le = int (1/ (f lx*dt) +0. 5) 
if (le.eq.2*(le/2) ) le = le+1 
le = le*3 + lb - 1 
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NFT=8192 
LEB=LE-LB+1 
PI=4 . 0*ATAN (1.0) 

DF=1 . 0/ (NFT*DT) 

PID=180. 0/PI 
do m = 1,37 

phlx = (xn-19)*10 
do n = 1,37 

ph2x = (n-19) *10 


sumax = 0.0 
sumf x = 0.0 
sumpx = 0.0 
sumaxe = 0.0 
sumf xe = 0.0 
sumpxe = 0.0 
amax =0.0 
f max = 0.0 

pmax = 0.0 , _ . 

CALL CREATE (xinit, A1X, FIX, PHlX,ph2x) 

Loop to create noise in data 
do k = 1,50 

DO I = 1,1000,2 

L vl = 2 . 0*ranl ( idum) - 1.0 

v2 = 2 . 0*ranl ( idum) - 1.0 
r = vl**2 + v2**2 
if (r.ge.l) go to 1 

fac = sqrt(- 2 . 0 *log(r)/r)*sd 
x(i) = vl*f ac + xinit ( i) 
x(i+i) = v2*f ac + xinit (i+1) 

CALL^WORK ( X , ampx , f rx , phx , AXer , PXer , FXer , A1X , FIX , PH1X) 

ax (k) = ampx 
fx(k) = frx 
px(k) = phx 
axe (k) = axer 
fxe(k) = fxer 
pxe(k) = pxer 
if (axer .gt. amax) 
if (fxer.gt. fmax) 
if (pxer. gt. pmax) 
sumax — sumax + ampx 
sumf x = sumfx + frx 
if (abs(phlx) .eq.180) then 

sumpx = sumpx + abs(phx) 


amax = axer 
fmax = fxer 
pmax = pxer 


else 

sumpx = sumpx + phx 
end if 

sumaxe = sumaxe + axer 
sumfxe = sumfxe + fxer 
sumpxe = sumpxe + pxer 
end do 

avgax = sumax/50.0 
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avgfx = sumf x/ 5 0.0 
avgpx = suxnpx / 50 . 0 
avgaxe = sumaxe/50 


avgfxe = 
avgpxe = 
sumax2 = 
sumf x2 = 
sumpx2 = 
sumaxe2 
sumfxe2 
sumpxe2 = 0.0 
do j = 1,50 

sumax2 = 
sumfx2 = 

if 


sumfxe/50 
sumpxe/50 
0.0 
0.0 
0.0 
= 0.0 
= 0.0 


0 

0 

,0 


- avgpx) **2 


suittax2 + (ax( j) - avgax)**2 
sumfx2 + (fx( j) - avgfx) **2 
(abs(phlx) .eq. 180) then 

sumpx2 - sumpx2 + (abs(px(j)) 

sis© 

sumpx2 = sumpx2 + (px(j) - avgpx) **2 
end if 

(axe ( 3 ) 

(fxe(j) 


= sumaxe2 
= sumfxe2 
= sumpxe2 


+ 

+ 

+ 


sumaxe2 
sumf xe2 
sumpxe2 
end do 

sdax = sqr t ( sumax2 / 4 9 . 0 ) 
sdfx = sqrt (sumf x2/ 49.0) 
sdpx = sqrt ( sumpx2 / 4 9 . 0 ) 
sdaxe = sqrt (sumaxe2/49 . 0) 
sdfxe = sqrt ( sumf xe2/ 4 9.0) 
sdpxe = sqrt ( sumpxe2 / 4 9 . 0 ) 
write ( 8 , * ) phlx , ph2x , amax 
write ( 9 , * ) phlx , ph2x , f max 
write(10, *)phlx,ph 2 x,pmax 
write (95,*) phlx , ph2x , avgaxe 
write (96 , * ) phlx , ph2x , avgfxe 
write (97 , * ) phlx , ph2x , avgpxe 
end do 
end do 
END 


(pxe( j) 


- avgaxe) **2 

- avgfxe) **2 

- avgpxe) **2 


SUBROUTINE CREATE (xt , A 1 X, FIX, PH1X, ph2x) 

REAL xt (3000) 

dt = 1.024 

iper =1000 

alx= 0. 02 

a0x=0 . 065 

a 2 x= 0 . 5 

pi = 4 . 0*ATAN (1.0) 
convrt = pi/ 180.0 
phi lx = phlx*convrt 
phi 2 x = ph 2 x*convrt 
TPIDT = 2 . 0 *PI*DT 
F2X = 0.03125 
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xtheta = tpidt*f lx 
xphi = tpidt*f2x 
DO I = 1 , IPER 

11 = rt(I) - AOX + A1X*C0S (xtheta*Il+PHilX) +A2X*COS (xphi*Il + PHi2X) 

END DO 
RETURN 
END 



IllCSSrrJSSSiErs" 

THE TIME INDEX WHERE THE MAXIMUM VALUE OCCURS. 

THIS IS BASED ON MODEL OF COS (WT+PHASE) . 

SUBROUTINE WORK ( ANG , amp , f req , phase , Aerr or , Perror , Ferror , A1 , FI ; PHI 1 ) 


PAR^METER D ^NDIM=3000 , NCDIM=8200 , NFT=8192 f 

DIMENSION AUX(NDIM) ,ANG(1) 


NPNT=7 ) 


REAL* 4 XFREQ (7) , PHIMAG (7 ) , PHREAL (7 ) 
COMPLEX AWO(NCDIM) 



COMMON / PAR 1 / LB , LE 
COMMON /PAR 2 /DT , PI , DF , PID 


T^ — FORCED TO BE ODD IN MAIN PROGRAM. 

HANN WINDOW ROUTINE USES ODD NUMBER OF POINTS TO TAP • 


LOAD INPUT DATA FROM ANG (I) INTO ARRAY AUX(J) . 
LB 1=1 “LB 
DO I=LB,LE 
IL=I+LB1 
AUX(IL)=ANG(I) 

END DO 


APPLY WINDOW FUNCTION TO TIME SEQUENCE 
CALL HANN (NTB1 , AUX , BIAS) 


MAKE COMPLEX NUMBER AWO(I) FROM REAL NUMBER AUX (I) 
HerH!£ginaot VALUE (AUX(I) IS the REAL VALUE). 


BY USING 


DO 1=1 / NTB1 

AWO ( I ) =CMPLX (AUX ( I ) , 0 . ) 
END DO 
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NOW PAD THE DATA STREAM WITH ZEROS OUT TO AWO(8192). 

DO I=NTB1+1 , NFT 

AWO(I) = CMPLX ( 0 . , 0 . ) 

END DO 


C SUBROUTINE FOUR1 DOES THE FOURIER TRANSFORM USING A FFT METHOD. 

C 

CALL FOUR1 ( AWO , NFT , 1 ) 

C LOOP TO FIND THE MAXIMUM MODULUS VALUE OF THE FOURIER TRANSFORM 
ampMAX=0 . 0 

istart — int ( (f 1 ~ 0 . 001) /df ) +1 
iend = int ( (f 1+0 . 001) /df ) +1 
do i = istart, iend 
fr = (i-1) *df 

if ( cabs ( awo ( i) ) .gt. ampmax) then 
ampmax = cabs(awo(i)) 
kf = i 
freq = fr 
end if 
end do 


C 

C 

C 

C 

C 

C 

C 

C 

C 


C 


C 


•RFATE THE 3 DATA SETS TO BE FITTED BY LEAST SQUARES POLYNOMIAL 
’OLYNOMIAL IS 2ND DEGREE AND 7 POINTS WILL BE USED IN CURVE FIT 
3 SETS ARE X 

MAGNITUDE OF TRANSFORM ( SQRT ( REAL* * 2 + IMAG**2)) 

REAL PART 
IMAGINARY PART 

; ENTER OF DATA SET IS THE FREQUENCY POINT WHERE MAX WAS FOUND. 


DO I = 1 , NPNT 

J = KF- ( (NPNT+1) /2 . 0) +1 
XFREQ(I) = CABS (AWO (J)) 
PHIMAG (I) = AIMAG(AWO(J) ) 
PHREAL(I) = REAL ( AWO ( J ) ) 
END DO 


DO CURVE FIT ON 
CALL TO LSCF WITH 
CURVE FITS AND 


THE MODULUS OF THE FOURIER TRANSFORM 
OPTION 1 DOES 2 THINGS. 

COMPUTES TRUE MAX FREQUENCY POINT. 


CALL LSCF (FQ_P0, ampMAX, XFREQ , 1) 

CALL LSCF (FQ_P0, PHASEI, PHIMAG, 2) 
CALL LSCF (FQ_P0, PHASER, PHREAL, 2) 
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non 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


FREQ = FREQ + DF * FQ_P0 

SCALING OF TRANSFORMED DATA IS PERFORMED TO GIVE OUTPUTS IN 
DEGS/SEC AND REPRESENT ACTUAL RATE DATA. 

SCALE = 4.0/ FLOAT (NTB1-1) 

AMP = SCALE * ampMAX 

PHASE = -AT AN 2 (PHASEI , PHASER) *pid 

THE FORWARD FOURIER ? SU ^E Y EXP^i*2*PI*F*T^ P ( THESE^WO* 

MANY FFT ROOTLES, ^CLUDING FOURl^UJ TRANSFO rm RESULT IN TWO 
DIFFERENT CONVENTIONS TH in the FIR st case, the SHIFT 

DIFFERENT FORMS FOR ^^ F transf0RM s AS G(F) , THEN G(T+T1) TRANSFORMS 

theorem states if G(T) transformS nd if g(t) TRANSFORM s as 

AS EXP ( i*2*PI*F*Tl) *G(F) • IN TH *9*PI*F*T1) *G (F) . SINCE OUR 

G(F), THEN G(T+T1) TRANSFORMS AS EXP- t PI p/^2*PI*F) ) , AND WE USE 

model is COS (2*PI*F*T + p) — SoupteR^TRANSFORM, we : expect our phase 
the first convention for the f °™ fe * r tr ^ nce the program uses the 

to BE 2*PI*F*P/(2*PI*F) - THE PHASE IS -P, SO TO 

correct'for^this^difference we must include another - sign: r (-P) = F * 

IS:$» a»°*Vo 100 

iSs(pnil) .‘S^P-o^U \IL (phase, -ahs (Phil) ) /360 . 0, *100 . 0 

return 

END 

SUBROUTINE HANN (LA, All, BIAS) 

SlI%=D C S=-F R OM INPUT SIGNAL 
PARAMETER (PI=3 . 14 15926) 

REAL All (LA) , HW(3000) 

ITM = (LA-1) / 2 
RM = FLOAT (ITM) 

DO IT= -ITM, ITM 

t — 1 + IT + ITM 

HW (I) =0 . 5* ( 1 • 0 + COS (PI* FLOAT ( IT) /RM ) ) 

A11(I)=AH(I) *HW(I) 

END DO 

CALL MEAN (LA, All, BIAS) 

BIAS = BIAS * 2.0 * LA/ (LA-1) 

D0 All(I) I = All (I) - BIAS * HW(I) 

END DO 
RETURN 
END 
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SUBROUTINE MEAN (LA, A2 2, SA) 


THIS ROUTINE COMPUTES 
MEAN IS NOT REMOVED, 


THE DC TERM OF THE DATA STREAM. 
BUT ONLY COMPUTED. 


REAL A22 (LA) 

SA = 0. 

DO 1=1, LA 

SA=SA+A22 (I) 
END DO 

SA=SA/FLOAT (LA) 

RETURN 

END 


SUBROUTINE FOUR1 (DATA, NN , I SIGN) 

C THIS ROUTINE DOES THE FOURIER TRANSFORM USING A FFT METHOD 

C REAL* 8 WR , WI , WPR , WP I , WTEMP , THETA 

DIMENSION DATA ( * ) 

N=2*NN 

J=1 

DO 11 1=1, N, 2 
IF(J.GT.I)THEN 
TEMPR=DATA(J) 

TEMPI=DATA( J+l) 

DATA ( J ) =DATA ( I ) 

DATA (J+l) =DATA(I+1) 

DATA ( I ) =TEMPR 
DATA(I+1)=TEMPI 
ENDIF 
M=N/2 

1 IF ( (M.GE.2) .AND. (J.GT.M) ) THEN 

J=J-M 
M=M/2 
GO TO 1 
ENDIF 
J=J+M 

11 CONTINUE 

MMAX=2 

2 IF (N.GT.MMAX) THEN 

ISTEP=2 *MMAX 

THETA=6 .283 18 5307 17 959 DO/ (ISIGN*MMAX) 

Wpr=-2 . D0*DSIN ( 0 . 5D0*THETA) **2 
WPI=DSIN (THETA) 

WR=1 . DO 
WI=0 . DO 

DO 13 M=1,MMAX, 2 
DO 12 I=M, N , I STEP 
J=I+MMAX 
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TEMPR=SNGL(WR) *DATA(J) -SNGL(WI) *DATA(J+1) 
TEMPI=SNGL(WR) *DATA(J+1) +SNGL(WI) *DATA( J) 
DATA(J) =DATA(I) -TEMPR 
DATA ( J+l) =DATA (1+1) -TEMPI 
DATA ( I ) =DATA ( I ) +TEMPR 

DATA(I+1)=DATA(I+1)+TEMPI 

CONTINUE 


WTEMP=WR 

WR=WR* WPR-WI * WP I +WR 
WI=WI *WPR+WTEMP *WP I+WI 


CONTINUE 

MMAX=ISTEP 


GO TO 2 
ENDIF 
RETURN 
END 



THIS SUBROUTINE DOES LEAST SQUARES CURVE FIT TO 7 

TOR A 2ND DEGREE POLYNOMIAL. THE DATA IS ASSUMED TO BE SAMPLED 

M INTEGRAL INTERVALS. ANY SCALING MUST BE DONE OUTSIDE 

THIS SUBROUTINE. THE 7 POINTS ARE : 

P = —3. —2, — 1 , 0, 1 , 2/ 3 
THE POLYNOMIAL IS F(P) = A + B * P + P * 

™IqSnc? C ?^IeSPONDING°tS Po'lS PO*DF (DF OF DATA STREAM) 
THIS DELTA IS REFERENCED TO MIDPOINT FREQUENCY OF 7 POIN . 
THE MAX VALUE IS F(PO) = A - (B*B) / (4*C) . 

SUBROUTINE LSCF (PO, FMAX, U_IN, IOPT) 



ON ENTRY: , 

U IN IS INPUT ORDINATE VALUES. (7 ) 

TAnm to OPTTON FOR 1 OF 2 THINGS 

1 : FIND PO WHERE MAX OCCURS PLUS COMPUTE MAX p0 
9 • COMPUTE VALUE OF POLYNOMIAL AT SPECIFIED FREQUNCY PO. 

IF IOPT=2 THEN PO IS FREQUENCY POINT TO EVALUATE POLYNOMIAL. 


PO IS VALUE OF P WHERE MAX PEAK OCCURS; 

THIS IS WRT CENTER POINT OF DATA. 
FMAX IS VALUE OF FUNCTION AT P=PO. 
******************************************* 


REAL* 4 U_IN ( * ) 

LOGICAL G FLAG 

********************************** 


****** 


FIRST STEP IS TO 
ALL COEFFICIENTS 
'A' IS -8, 12, 
'B' IS -9, -6, 
•C IS 5, 0, -3 


DO LEAST SQUARES. 

HAVE BEEN PRE-COMPUTED . 
24, 28, 24, 12, -8 DIVIDED 
-3, 0, 3, 6, 9 DIVIDED 

, -4, -3, 0, 5 DIVIDED 


BY 

BY 

BY 


84 

84 

84 
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US17 = U_IN ( 1) + U_IN (7) 

US35 = U_IN ( 3 ) + U_IN ( 5 ) 

A = C0F1 

COFl = -8 . *US17 + 12 . * (U_IN ( 2 ) +U_IN ( 6 ) ) + 

# 24 . *US3 5 + 28 . *U_IN (4) 

B = COF2 

COF2 = 9.*(“U IN ( 1) +U_IN (7 ) ) + 

# 6 . * ( -U_IN ( 2 ) +U_IN ( 6 ) ) + 

f 3 . * ( -U_IN ( 3 ) +U_IN ( 5 ) ) 

0Qp3 

COF3 = 5 . *US17 -3 . *US35 - 4.*U_IN(4) 

IF ( ABS ( COF3 ) . LT . 1 . 0E-08 ) THEN 

PRINT* ' ********* WARNING ********* 

PRINT* ' CONSTANT VALUE EQUALS ZERO - CANNOT COMPUTE A 
PRINT* ,' FREQUENCY VALUE.' 

RETURN 

ENDIF 

DID NOT DIVIDE BY 84 YET. DO SO FOR MAX PART BUT NOT PO. 

COMPUTE PO, VALUE OF F(P) WHERE A+B*P+ C * P *P = 0 
COMPUTE FUNCTION AT PO; A+B*P0+C*P0*P0 - A-B 2/4 C 
IF ( IOPT .EQ. 1) THEN 
PO=-0.5*COF2/COF3 

FMAX = (COFl + 0.5 * PO * COF2) / 84. 


ELSE 

FMAX = (COFl + 
END IF 
RETURN 
END 


P0*( COF2 + PO*COF3 ) ) /84 


ial = 7141, 
ia2 = 8121, 
ia3 = 4561, 


icl = 54773, rml = 1 
ic2 = 28411, rm2 = 1 
ic3 = 51349) 


function ranl(idum) 
dimension r(97) 
parameter (ml = 259200, 

parameter (m2 = 134456, 

parameter (m3 = 243000, 

data iff /0 / 

if ( idura. It . 0 . or . if f • e<3« 0) then 

iff = 1 

ixl = mod (icl - idum,ml) 

ixl = mod(ial*ixl + icl, ml) 

ix2 = mod (ixl, m2) 

ixl = mod ( ial*ixl + icl, ml) 

ix3 = mod (ixl, m3) 

do j = 1,97 . 

ixl = mod(ial*ixl + icl,ml) 

ix2 = mod ( ia2*ix2 + ic2,m2) 

r ( j ) = (float (ixl) + f loat (ix2) *rm2) *rml 

end do 
idum = 1 
end if 


MAXIMUM' 


.0/ml) 
. 0/m2) 
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ixi = mod ( is 1 * ixl + id t itil) 
ix2 = mod( ia2*ix2 + ic2 # iu2) 
ix3 — mod ( ia3 * ix3 + ic3,m3) 
j » 1 + (97*ix3) /m3 
if ( j .gt.97 .or. j . It. 1) pause 

rtj)" = Afloat (ixl) + float (ix 2 ) *rm2) *rml 

return 

end 


C Program LIBRATION . FOR SSh^mIaH^FODRI^ '^ and^IcF subroutines 

= f rom^the^UNOMSC . FOR^prograra . ‘"libRATION^FOR syst4matically runs through the 
C phase pairings for a given frequency. 


REAL X(4000) 

COMMON / PARI / LB , LE , NFT , NPNT , NDEG 
COMMON/ PAR2 / DT , PI , DF , PID 
dt = 1.024 
sd = 2 . 8e-03 


lb=l 


npnt = 7 

ndeg = 3 

fix = 0.0019 

read* , numcycle , iseed 

le = int (1/ (f lx*dt) +0.5) 

if (le.eq. 2* (le/2) ) le = le+1 

le = le*numcycle + lb - 1 


NFT=8 192 
LEB=LE-LB+1 
PI=4 . 0*ATAN (1.0) 

DF=1.0/(NFT*DT) 

PID=180. 0/PI 
do m = 1/37 

phlx = (m-19) *10 
do n = 1/37 

ph2x = (n-19) *10 

CALL CREATE (X,A1X, FIX, PHlX,ph2x) 

CALL WORK ( X , ampx , f rx , phx , aerror , 
write ( 11 , * ) phlx , ph2x, aerror 
write (12 , *) phlx,ph2x, ferror 
wr ite ( 13 , * ) phlx , ph2x , perror 
end do 
end do 


ferror , perror ,A1X, FIX, PH1X) 


END 
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SUBROUTINE CREATE (X , A1X , FIX , PHI IX, phi2x) 

REAL X(4000) 
tlngth=20000 . 0 

dt = 1.024 
iper=3000 
alx=0 . 0034 
a0x=0 . 065 
a2x=0.05 
alibx = 0.0046 
PI = 4 . 0*ATAN (1.0) 
convrt = pi/ 180.0 
phlx = philx*convrt 
ph2x = phi 2 x*convrt 
TPIDT = 2 . 0*PI*DT 
fix = 0.0019 
F2X = 0.089 
flibx = 1/2713.0 
xtheta = tpidt*f lx 
xphi = tpidt*f2x 
xlibr = tpidt*flibx 
DO I = 1 , IPER 

II = 1-1 

x(I) S ’- 1 A0X + A1X*C0S (xtheta*Il+PHlX) +A2X*COS (xphi*Il) + 
$ alibx*cos (xlibr*il+ph2x) 

END DO 
RETURN 
END 


C 

C 

c 

c 

c 

c 

c 


SUBROUTINE WORK CALCULATES THE AMPLITUDE, PHASE, AND FREQUENCY 
O? THE DATA THE FOURIER TRANSFORM SUBROUTINE FOUR1 IS 
^ATTFD BY SUBROUTINE WORK. WORK RETURNS TO THE MAIN PROGRAM 
THE** VALUES OF THE AMPLITUDE, PHASE, AND FREQUENCY AS WELL AS 
THE TIME INDEX WHERE THE MAXIMUM VALUE OCCURS. 

THIS IS BASED ON MODEL OF COS (WT+PHASE) . 

;UBROUTINE WORK ( ANG , amp , f r eq , phase , Aerror , Perr or , Ferr or , A1 , F 1 , PHI 1 ) 


INTEGER NDIM, NCDIM 

PARAMETER (NDIM=4000, NCDIM=8200, 
DIMENSION AUX (NDIM) , ANG (1) 


NFT=8 192 , NPNT=7 ) 


REAL* 4 XFREQ(7) , PHIMAG ( 7 ) , PHREAL ( 7 ) 
COMPLEX AWO (NCDIM) 


COMMON / PARI / LB , LE 

COMMON/ PAR2 / DT , PI , DF , PID 


C 

C 


nBl IS FORCED TO BE ODD IN MAIN PROGRAM. 

yjN WINDOW ROUTINE USES ODD NUMBER OF POINTS TO TAPER. 
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LOAD INPUT DATA FROM ANG(I) INTO ARRAY AUX(J) . 
LB 1=1~ LB 
DO I=LB,LE 
IL=I+LB1 
AUX(IL)=ANG(I) 

END DO 


APPLY WINDOW FUNCTION TO TIME SEQUENCE 
CALL HANN (NTB1 , AUX , BIAS) 



make complex number AWO(I) from keal NUMBER AUX (I) BY using 
a ZERO IMAGINARY VALUE (AUX(I) IS THE REAL VALUE). 


DO 1=1 , NTB1 

AWO ( I ) =CMPLX (AUX (I) ,0.) 
END DO 


C NOW PAD THE DATA STREAM WITH ZEROS OUT TO AWO(8192) . 

C 

DO I=NTB1+1 , NFT 

AWO (I) = CMPLX ( 0 . ,0. ) 

END DO 

C SUBROUTINE FOUR1 DOES THE FOURIER TRANSFORM USING A FFT METHOD 

C 

CALL FOURl(AWO,NFT,l) 

C LOOP TO FIND THE MAXIMUM MODULUS VALUE OF THE FOURIER TRANSFORM 
ampMAX=0 . 0 

istart = int( (f 1 - 0 . 001) /df ) +1 
iend = int ( (f 1+0. 001) /df ) +1 
do i = istart, iend 
fr = (i-1) *df 

if (cabs (awo (i) ) *gt .ampinax) then 
ampmax = cabs (awo (i)) 
kf = i 
freq = fr 
end if 
end do 
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c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

_c 

c 

c 

c 

c 

c 

c 

-c 

c 

c 

c 


CREATE THE 3 DATA SETS TO BE FITTED BY LEAST SQUARES POLYNOMIAL. 
POLYNOMIAL IS 2ND DEGREE AND 7 POINTS WILL BE USED IN CURVE FIT. 

3 SETS ARE • 

MAGNITUDE OF TRANSFORM ( SQRT (REAL**2 + IMAG**2) ) 

REAL PART 
IMAGINARY PART 

CENTER OF DATA SET IS THE FREQUENCY POINT WHERE MAX WAS FOUND. 


DO I = 1 , NPNT 

J = KF-( (NPNT+1) /2.0)+I 
XFREQ(I) = CABS ( AWO ( J) ) 
PHIMAG ( I ) = AIMAG (AWO ( J) ) 
PHREAL ( I ) = REAL (AWO ( J) ) 

END DO 


DO CURVE FIT ON THE MODULUS OF 
CALL TO LSCF WITH OPTION 1 DOES 
CURVE FITS AND COMPUTES TRUE 


THE FOURIER TRANSFORM 
2 THINGS. 

MAX FREQUENCY POINT. 


CALL LSCF (FQ_PO, ampMAX, XFREQ , 1) 


CALL LSCF (FQ_PO, PHASEI, PHIMAG, 2) 
CALL LSCF (FQ_PO , PHASER, PHREAL, 2) 
FREQ = FREQ + DF * FQ_PO 


SCALING OF TRANSFORMED DATA IS PERFORMED TO GIVE OUTPUTS IN 
DEGS/SEC AND REPRESENT ACTUAL RATE DATA. 

SCALE = 4 . 0/ FLOAT (NTB1-1) 

AMP = SCALE * ampMAX 

PHASE = -ATAN2 (PHASEI , PHASER) *pid 


THE FORWARD FOURIER TRANSFORM IS USUALLY DEFINED WITH EXP ( -1*PI*F*T) . 
MANY^FFT ROUTINES, INCLUDING FOUR1, USE EXP (+i*2*PI*F*T) . THESE TWO 
DIFFERENT CONVENTIONS FOR THE FORWARD FOURIER TRANSFORM RESULT IN TW° 
nTEEFRENT FORMS FOR THE SHIFT THEOREM. IN THE FIRST CASE, THE SHIFT 
THEOREM STATES THAT IF G (T) TRANSFORMS AS G(F) , THEN G(T+T1) TRANSFORMS 
AS EXWi*2*PI*F*Tl *G(F) IN THE SECOND CASE, IF G(T) TRANSFORMS AS 
cm TOEN G(T+T1) TRANSFORMS AS EXP(-i*2*PI*F*Tl)*G(F) . SINCE OUR 
Sonk IS COS 2*PI.F*T + P) = COS (2*PI*F*(T + P/(2*PI*F>), AND WE USE 
THE FIRST CONVENTION FOR THE FOURIER TRANSFORM, WE EXPECT OUR P |^ SE 
TO BE 2*k*F^/ (2.PI.F) = P. HOWEVER, SINCE THE PROGRAM USES THE 
SECOND CONVENTION FOR THE FOURIER TRANSFORM, THE PHASE IS -P , SO TO 
CORRECT FOR THIS DIFFERENCE WE MUST INCLUDE ANOTHER - SIGN. -( P) - P- 


FERROR 

AERROR 

perror 

RETURN 

END 


(ABS (FREQ-F1) /0 . 005) * 100 
(ABS (AMP-A1) /Al) * 100 

(abs (abs (phase) -abs (phil) ) / 360. 0) *100.0 
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c 

c 

c 



SUBROUTINE HANN (LA, All, BIAS) 


TAPER IS RAISED COSINE CURVE. 

MEAN IS COMPUTED AND REMOVED FROM INPUT SIGNAL 
PARAMETER (PI=3 . 1415926) 

REAL All (LA) , HW(3000) 

ITM = (LA-1) / 2 
RM = FLOAT (ITM) 

DO IT= -ITM, ITM 

T = 1 + IT + ITM . . 

HW ( I) = 0 . 5* ( 1 • 0 + COS (PI*FLOAT ( IT) /RM ) ) 

All (I) “All (I) *HW (I) 

END DO 

™ .... 


CALL MEAN (LA, All, BIAS) 

BIAS = BIAS * 2.0 * LA/ (LA-1) 

DO I = 1/LA , , 

All (I) = All ( I ) “ BIAS * HW (I) 

END DO 
RETURN 
END 


SUBROUTINE MEAN (LA, A2 2, SA) 


THIS ROUTINE COMPUTES 
MEAN IS NOT REMOVED, 


THE DC TERM OF THE DATA STREAM. 
BUT ONLY COMPUTED. 


REAL A22 (LA) 

SA = 0. 

DO 1=1, LA 

SA=SA+A22 (I) 
END DO 

SA=SA/ FLOAT (LA) 

RETURN 

END 


SUBROUTINE FOUR1 (DATA, NN, ISIGN) 

THIS ROUTINE DOES THE FOURIER TRANSFORM USING A FFT METHOD 

REAL* 8 WR, WI,WPR, WPI ,WTEMP, THETA 
DIMENSION DATA(*) 

N=2*NN 

J=1 

DO 11 1=1, N, 2 
IF(J.GT.I)THEN 
TEMPR=DATA ( J ) 
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TEMPI=DATA (J+l) 

DATA ( J ) =D ATA ( I ) 

DATA (J+l) =DATA(I+1) 

DATA ( I ) =TEMPR 
DATA(I+1)=TEMPI 
ENDIF 
M=N/2 

I IF ( (M.GE.2) .AND. (J.GT.M) ) THEN 

J=J-M 
M=M/2 
GO TO 1 
ENDIF 
J=J+M 

II CONTINUE 
MMAX=2 

2 IF (N.GT.MMAX) THEN 

ISTEP=2*MMAX 

THETA=6 .28318530717959D0/ (ISIGN*MMAX) 
WPR=-2.D0*DSIN(0.5D0*THETA) **2 
WPI=DSIN (THETA) 

WR=1 • DO 
WI=0.D0 

DO 13 M=1 , MMAX , 2 
DO 12 I=M, N , I STEP 

J=I+MMAX , , , % 

TEMPR=SNGL(WR) *DATA(J) -SNGL(WI) *DATA( J+l) 

TEMPI=SNGL (WR) *DATA ( J+l) +SNGL (WI) *DATA ( J) 

DATA ( J) =DATA ( I ) -TEMPR 

DATA (J+l) =DATA(I+1) -TEMPI 

DATA ( I ) =DATA ( I ) +TEMPR 

DATA ( 1+1 ) =DATA ( 1+1 ) +TEMPI 

12 CONTINUE 
WTEMP=WR 

WR=WR*WPR-WI *WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 

13 CONTINUE 
MMAX=ISTEP 

GO TO 2 
ENDIF 
RETURN 
END 


C 

C 

C 

C 

C 

C 

C 

C 

C 


HIS SUBROUTINE DOES LEAST SQUARES CURVE FIT TO 7 POINTS 
OR A 2ND DEGREE POLYNOMIAL. THE DATA IS ASSUMED TO BE SAMPLED 
T INTEGRAL INTERVALS. ANY SCALING MUST BE DONE OUTSIDE 
HIS SUBROUTINE. THE 7 POINTS ARE : 

P = -3 1 “2, “* 1 / 0 , 1/ 2, 3 

HE POLYNOMIAL IS F(P) = A + B*P + C*P*P. 

FREQOTNC? C C^IeSPONDING 0 TO P0^IS*P0*DF (DF OF DATA STREAM) 
?HIS DELTA IS REFERENCED TO MIDPOINT FREQUENCY OF 7 POINTS. 
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THE MAX VALUE IS F(PO) = A - (B*B)/(4*C). 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

C 

C 

C 

C 


c 


c 


c 


SUBROUTINE LSCF (PO, FMAX, U_IN, IOPT) 


ON ENTRY : 

U IN IS INPUT ORDINATE VALUES. (7 ) 

IOPT IS OPTION FOR 1 OF 2 THINGS 

1 : FIND PO WHERE MAX OCCURS PLUS COMPUTE MAX VALUE. 

2 : COMPUTE VALUE OF POLYNOMIAL AT SPECIFIED FREQUNCY PO. 
IF IOPT=2, THEN PO IS FREQUENCY POINT TO EVALUATE POLYNOMIAL. 


ON EXIT 

PO IS VALUE OF P WHERE MAX PEAK OCCURS; 

THIS IS WRT CENTER POINT OF DATA. 

FMAX IS VALUE OF FUNCTION AT P=PO . 
******************************************** 

REAL* 4 U_IN ( * ) 

LOGICAL G FLAG 

******************************************** 

FIRST STEP IS TO DO LEAST SQUARES. 

ALL COEFFICIENTS HAVE BEEN PRE-COMPUTED. 

'A' IS -8, 12, 24, 28, 24, 12, -8 DIVIDED BY 84 

*B' IS -9, -6, -3, 0, 3, 6, 9 DIVIDED BY 84 

'C' IS 5, 0, -3, -4, -3, 0, 5 DIVIDED BY 84 


US17 = U_IN (1) + U_IN (7) 

US35 = U_IN ( 3 ) + U_IN ( 5 ) 

A = COF1 

COF1 = -8 . *US17 + 12 . * (U_IN (2) +U_IN (6) ) + 

/ 24 . *US35 + 28 . *U_IN (4) 

B = COF2 

COF2 = 9 . * ( -U IN ( 1 ) +U_IN ( 7 ) ) + 

# 6 . * ( -U_IN ( 2 ) +U_IN ( 6 ) ) + 

# 3 . * ( -U_IN ( 3 ) +U_IN (5) ) 

C= COF3 

COF3 = 5 . *US17 -3 . *US35 - 4.*U_IN(4) 

IF (ABS(COF3) .LT.1.0E-08) THEN 

PRINT* ********** WARNING ********** 

PRINT*, 'CONSTANT VALUE EQUALS ZERO - CANNOT COMPUTE A MAXIMUM' 
PRINT* \ ' FREQUENCY VALUE.' 

RETURN 

ENDIF 


DID NOT DIVIDE BY 84 YET. DO SO FOR MAX PART BUT NOT PO . 

COMPUTE PO, VALUE OF F(P) WHERE A+B*P+C*P*P = 0 
COMPUTE FUNCTION AT PO; A+B*P0+C*P0*P0 = A-B**2/4C 
IF (IOPT .EQ. 1) THEN 
P0=-0 . 5*COF2/COF3 

FMAX = (COF1 + 0.5 * PO * COF2) /84. 


ELSE 
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FMAX = (C0F1 + P0*( COF2 + PO*COF3 ) )/84. 

END IF 
RETURN 
END 


c Program LIB_NOISE . FOR uses the signal and noise generation algorithms from 
C^RELIBR.FOR program and the WORK, HANN, MEAN, FOUR1, and LSCF subroutines 

from TT _ MnTqp FOR systematically runs through the phase 

c the UNOMSC . FOR program. LIB_NOISE . FOR syscema * 

C pairings for a given frequency. 

REAL X ( 3000) ,xinit(3000) , ax ( 100 ) , fx ( 100) ,px(100) ,axe(100) .fxe<100) , 
i pxe(100) 

EXTERNAL rani m 

COMMON /PARI /LB , LE , NFT , NPNT , NDEG 
COMMON/ PAR2 / DT , PI , DF, PID 

dt = 1.024 

sd = 2 . 8e — 03 

read* , numcycle , iseed 

lb=l 

npnt = 7 
ndeg = 3 
fix = 0.0019 

le = int(l/ (flx*dt) +0.5) 
if (le. eq. 2* (le/2) ) le = le+1 
le = le*numcycle + lb - 1 
NFT=8192 
LEB=LE-LB+1 
PI=4 . 0*ATAN (1.0) 

DF=1.0/(NFT*DT) 

PID=180. 0/PI 
do m = 1 r 37 

phlx = (m-19 ) *10 
do n = 1,37 

ph2x = (n-19) *10 

sumax = 0.0 
sumf x = 0.0 
sumpx = 0.0 
sumaxe = 0.0 
sumfxe = 0.0 
sumpxe = 0.0 
amax = 0.0 
fmax = 0.0 

CALL CREATE (xinit , A1X , FIX, PH1X , ph2x) 
c Loop to create noise in data 
do k = 1,50 
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DO I = 1,1000,2 

vi = 2 . 0*ranl ( idum) - 1.0 

v2 = 2 . 0*ranl (idum) - 1.0 

r = vl**2 + v2**2 

if (r.ge.l) go to 1 

fac = sqrt (-2 . 0*log (r) /r) *sd 

x(i) = vl*fac + xinit(i) 

x(i+l) = v2*f ac + xinit ( i+1) 

CALL WORK ( X , ampx , f rx , phx , AXer , PXer , FXer , A1X , FIX , PH1X) 

ax(k) = ampx 
fx(k) = frx 
px(k) = phx 
axe (k) = axer 
fxe(k) = fxer 
pxe(k) = pxer 
if (axer .gt.amax) 
if (fxer. gt.f max) 
if (pxer.gt.pmax) 
sumax = sumax + ampx 
sumfx = sumfx + frx 
if (abs (phlx) . eq. 180) then 

sumpx = sumpx + abs (phx) 

else 

sumpx = sumpx + phx 
end if 

sumaxe = sumaxe + axer 
sumfxe = sumfxe + fxer 
sumpxe = sumpxe + pxer 
nd do 

vgax = sumax/ 5 0.0 
vgfx = sumfx/ 50.0 
vgpx = sumpx/ 5 0.0 

sumaxe/ 50.0 
sumfxe/50 . 0 
sumpxe/ 50 . 0 
0.0 


amax = axer 
fmax = fxer 
pmax = pxer 


vgaxe 

vgfxe 

vgpxe 

umax2 

umfx2 

umpx2 


0.0 
0.0 

^ 

amaxe2 = 0.0 
imfxe2 = 0.0 
ampxe2 = 0.0 
3 j = 1/50 

sumax2 = sumax2 


sumaxe — oumuAi + (ax( j) avgax) *2 
sumf x2 = sumf x2 + (fx(j) - a vgfx) **2 
if (abs (phlx) . eq. 180) then , 

sumpx2 = sumpx2 + (abs(px(j)) avgpx) 2 

6 XS6 

sumpx2 = sumpx2 + (px(j) - avgpx) **2 

sumaxe2 = sumaxe2 + (axe(j) - avgaxe) **2 
sumfxe2 = sumfxe2 + (fxe(3) - avgfxe)**2 
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sumpxe2 = sumpxe2 + (pxe(j) - avgpxe)**2 
end do 

sdax = sqrt ( sumax2 / 4 9 . 0 ) 
sdfx = sqrt ( sumf x2 / 4 9 . 0 ) 
sdpx = sqrt ( sumpx2 / 4 9 . 0 ) 
sdaxe = sqrt (sumaxe2/49 . 0) 
sdfxe = sqrt(sumfxe2/49.0) 
sdpxe = sqrt (sumpxe2/49 . 0) 
write ( 8 , * ) phlx , ph2x , amax 
writs ( 9 i * ) phlx , ph2x , f max 

write(l° / *)P hlx ' ph2x,pinaX 

write ( 95 , * ) phlx , ph2x , avgaxe 

write(96, *)phlx,ph2x,avgfxe 

write ( 97 , * ) phlx , ph2x , avgpxe 
end do 
end do 
END 

SUBROUTINE CREATE (x,AlX,FlX,PHHX,phi2x) 

REAL x ( 3 000 ) f y (3000) 
tlngth=20000 . 0 
dt = 1.024 
iper=3000 
alx=0 . 0034 
a0x=0 . 065 
a2x=0 . 05 
alibx = 0 . 0046 
pi = 4 . 0*ATAN ( 1 . 0) 
convrt = pi/ 180.0 
phlx = philx*convrt 
ph2x = phi 2 x*convrt 
TPIDT = 2.0*PI*DT 
fix = 0.0019 
F2X = 0.089 
flibx = 1/2713.0 
xtheta = tpidt*f lx 
xphi = tpidt*f2x 
xlibr = tpidt*flibx 
DO I = 1 , IPER 
II = 1-1 

l™ = X A0X + A1X*C0S (xtheta*Il+PHlX) +A2X*COS (xphi*Il) + 

$ alibx*cos (xlibr*il+ph2x) 

END DO 
RETURN 
END 


25 



_ . 7ADV raTCULATES the amplitude, phase, and frequency 

&*»* 

TOE L VALUES S of?HE MPLITUDE, PHASE, AND FREQUENCY AS WELL AS 

the time index where the maximum value occurs, 
this is based on MODEL OF COS (WT+PHASE) . 

SUBROUTINE WORK (ANG, amp, freq, phase, Aerror, Perror, Terror ,A1, FI, PHI1) 

S?ER D ™DlS. NCDIM=8200 , NFT-8I92. NPNT=7, 

DIMENSION AUX (NDIM) , ANG ( 1) 

REAL* 4 XFREQ ( 7 ) , PHIMAG ( 7 ) , PHREAL ( 7 ) 

COMPLEX AWO(NCDIM) 

COMMON / PARI / LB , LE 
COMMON / PAR2 /DT , PI , DF , PID 

EsffK " ~ 

LOAD INPUT DATA FROM ANG(I) INTO ARRAY AUX(J). 

LB 1=1 -LB 
DO I=LB,LE 
IL=I+LB1 
AUX(IL)=ANG(I) 

END DO 


apply window function to time sequence 

CALL HANN (NTB1 , AUX , BIAS) 

DO 1=1 , NTB1 

AWO(I)=CMPLX(AUX(I) , 0. ) 

END DO 

NOW PAD THE DATA STREAM WITH ZEROS OUT TO AWO(8192) . 

DO I=NTB1+1 , NFT 

AWO(I) = CMPLX ( 0 . ,0. ) 

END DO 
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non o oooo 


C SUBROUTINE F0UR1 DOES THE FOURIER TRANSFORM USING A FFT METHOD 

C 

CALL FOUR1 (AWO, NFT , 1) 

c LOOP TO FIND THE MAXIMUM MODULUS VALUE OF THE FOURIER TRANSFORM 
ampMAX=0 . 0 

istart = int ( (f 1-0 . 001) /df ) +1 
iend = int ( (f 1+0. 001) /df ) +1 
do i = istart, iend 
fr = (i-1) *df 

if ( cabs ( awo ( i ) ) . gt . ampmax ) then 
ampmax = cabs ( awo (i)) 
kf = i 
freq = fr 
end if 
end do 


C 

C 

C 

c 

c 

c 

c 

c 

c 


c 


POLYNOMIAL S5LTJZ 

3 S MAGNITUDE OF TRANSFORM (SQRT (REAL**2 + IMAG**2)) 

REAL PART 
IMAGINARY PART 

CENTER OF DATA SET IS THE FREQUENCY POINT WHERE MAX WAS FOUND. 


DO I = 1,NPNT 

J = KF- ( (NPNT+1) /2 . 0) +1 
XFREQ(I) = CABS (AWO (J)) 
PHIMAG(I) = AIMAG (AWO ( J) ) 
PHREAL ( I ) = REAL ( AWO ( J ) ) 
END DO 


)0 CURVE FIT ON THE MODULUS OF THE FOURIER TRANSFORM 
'an TO LSCF WITH OPTION 1 DOES 2 THINGS. 

CU^?E FITS AND COMPUTES TRUE MAX FREOUENCY POINT. 


CALL LSCF (FQ_P0 , ampMAX, XFREQ , 1) 

CALL LSCF (FQ_P0 , PHASEI, PHIMAG, 2) 

CALL LSCF (FQ_P0, PHASER, PHREAL, 2) 

FREQ = FREQ + DF * FQ_P0 

SCALING OF TRANSFORMED DATA IS PERFORMED TO GIVE OUTPUTS IN 
DEGS/SEC AND REPRESENT ACTUAL RATE DATA. 

SCALE = 4 . 0/FLOAT (NTB1-1) 

AMP = SCALE * ampMAX 

PHASE = -ATAN2 (PHASEI, PHASER) *pid 
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noon n n o 


_ c 
c 
c 
c 

~ c 
c 
c 

- c 
c 
c 

_ c 
c 
c 


cru* rnuwARn FOURIER TRANSFORM IS USUALLY DEFINED WITH EXP ( -i*PI*F*T) . 

SSS™*353£ioS^^ 

SF Lk%e y.s'issst j&.’ss s^rass- 

MODEi if COS ^PlUffp)- COS (2*PI*F* (T + P/(2*PI*F)), AND WE USE 
raE FlIsT CONTENTION FOR THE FOURIER TRANSFORM, WE EXPECT OUR PHASE 
TO BE 2*PI*F*P/ (2*PI*F) = P. HOWEVER, SINCE THE PROGRAM USES THE 
SECOND CONVENTION FOR THE FOURIER TRANSFORM, THE PHASE IS -P, SO TO 
CORRECT FOR THIS DIFFERENCE WE MUST INCLUDE ANOTHER - SIGN. -( P) p * 


FERROR = 
AERROR = 
perror = 
RETURN 
END 


(ABS (FREQ-F1) /0.005) * 100 

(ABS (AMP-A1) /Al) * 100 „ 

(abs(abs (phase) -abs(phil) ) / 3 60.0) *100.0 


SUBROUTINE HANN (LA, All, BIAS) 


TAPER IS RAISED COSINE CURVE. „ T ™ T 

MEAN IS COMPUTED AND REMOVED FROM INPUT SIGNAL 
PARAMETER (PI=3 . 1415926) 

REAL All (LA) , HW(3000) 

ITM = (LA-1) /2 
RM = FLOAT (ITM) 

DO IT= -ITM, ITM 

I = 1 + IT + ITM 

HW(I) =0 . 5* ( 1 . 0 + COS (PI* FLOAT ( IT) /RM ) ) 

A11(I)=A11(I) *HW(I) 


END DO 

COMPUTE MEAN OF TAPERED SIGNAL 
TRUE MEAN IS TWICE COMPUTED VALUE 
REDUCES VALUE BY FACTOR OF 2 . (I. 


BECAUSE HANN WINDOW 
E. MEAN OF WINDOW IS 


0.5) 


CALL MEAN (LA, All, BIAS) 

BIAS = BIAS * 2.0 * LA/ (LA-1) 

DO I = 1, LA 

All ( I ) = All ( I ) - BIAS * HW (I) 
END DO 
RETURN 
END 
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o o o o n n o 


SUBROUTINE MEAN (LA, A2 2 , SA) 

THIS ROUTINE COMPUTES THE DC TERM OF THE DATA STREAM. 
MEAN IS NOT REMOVED, BUT ONLY COMPUTED. 

REAL A22 (LA) 

SA = 0. 

DO 1=1, LA 

SA=SA+A2 2 (I) 

END DO 

SA=SA/ FLOAT (LA) 

RETURN 

END 


SUBROUTINE FOUR1 (DATA, NN , ISIGN) 

THIS ROUTINE DOES THE FOURIER TRANSFORM USING A FFT METHOD. 

REAL* 8 WR , WI , WPR , WPI , WTEMP , THETA 
DIMENSION DATA ( * ) 

N=2*NN 
J=1 

DO 11 1=1, N, 2 
IF ( J . GT . I) THEN 
TEMPR=DATA ( J ) 

TEMPI=DATA ( J+l) 

DATA( J) =DATA(I) 

DATA (J+l) =DATA ( 1+1) 

DATA ( I ) =TEMPR 
DATA ( I + 1 ) =TEMPI 
END IF 
M=N/2 

IF ( (M.GE.2) .AND. (J.GT.M) ) THEN 
J=J-M 
M=M/2 
GO TO 1 
ENDIF 
J=J+M 
11 CONTINUE 

MMAX=2 

2 IF (N.GT.MMAX) THEN 

ISTEP=2*MMAX 

THETA=6 .28318530717959D0/ (ISIGN*MMAX) 

WPR=-2 . D0*DSIN ( 0 . 5D0*THETA) **2 
WPI=DSIN (THETA) 

WR=1. DO 
WI=0 . DO 

DO 13 M=1,MMAX, 2 
DO 12 I=M, N , I STEP 

J=I+MMAX , , „ . 

TEMPR=SNGL (WR) * DATA ( J ) -SNGL(WI) *DATA(J+1) 
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13 


TEMPI=SNGL (WR) * DATA (J+l) +SNGL (WI) *DATA ( J) 
DATA ( J) =DATA ( I ) -TEMPR 
DATA ( J+l )=DATA(I+1) -TEMPI 
DATA ( I ) =DATA ( I ) +TEMPR 
DATA ( 1+1 ) =DATA ( 1+1 ) +TEMPI 
CONTINUE 
WTEMP=WR 

WR=WR*WPR-WI *WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 
CONTINUE 
MMAX=ISTEP 
GO TO 2 
END IF 
RETURN 
END 


C THIS SUBROUTINE DOES LEAST SQUARES CURVE FIT TO 7 POINTS 

C FOR A 2ND DEGREE POLYNOMIAL. THE DATA IS ASSUMED TO BE SAMPLED 

C AT INTEGRAL INTERVALS. ANY SCALING MUST BE DONE OUTSIDE 

C THIS SUBROUTINE. THE 7 POINTS ARE : 

C P = -3, -2, -1, 0, 1, 2, 3 

C THE POLYNOMIAL IS F(P) = A + B*P + C*P*P. 

C THE MAX OCCURS AT P = P0 = -B/(2*C). 

C FREQUENCY CORRESPONDING TO P0 IS P0*DF (DF OF DATA STREAM) 

C THIS DELTA IS REFERENCED TO MIDPOINT FREQUENCY OF 7 POINTS. 

C THE MAX VALUE IS F(P0) = A - (B*B)/(4*C). 

C 

SUBROUTINE LSCF (P0, FMAX, U_IN, IOPT) 


C 

c ON ENTRY : 

C U_IN IS INPUT ORDINATE VALUES. (7 ) 

C IOPT IS OPTION FOR 1 OF 2 THINGS 

C 1 : FIND P0 WHERE MAX OCCURS PLUS COMPUTE MAX VALUE. 

C 2 : COMPUTE VALUE OF POLYNOMIAL AT SPECIFIED FREQUNCY P0. 

C IF IOPT=2, THEN P0 IS FREQUENCY POINT TO EVALUATE POLYNOMIAL. 

C 

C ON EXIT 

C P0 IS VALUE OF P WHERE MAX PEAK OCCURS; 

C THIS IS WRT CENTER POINT OF DATA. 

C FMAX IS VALUE OF FUNCTION AT P=P0. 

Q ******************************************** 

REAL* 4 U_IN ( * ) 

LOGICAL G_FLAG 

C ******************************************** 

C FIRST STEP IS TO DO LEAST SQUARES. 

C ALL COEFFICIENTS HAVE BEEN PRE-COMPUTED . 

C 'A' IS -8, 12, 24, 28, 24, 12, -8 DIVIDED BY 84 

C 'B' IS -9, -6, -3, 0, 3, 6, 9 DIVIDED BY 84 

C 'C' IS 5, 0, -3, -4, -3, 0, 5 DIVIDED BY 84 

C 

US17 = U IN (1) + U_IN (7) 
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c 

c 

c 

c 

c 


US35 = U_IN ( 3 ) + U_IN ( 5 ) 

A = C0F1 

COFl = -8 . *US17 + 12 . * (U IN ( 2 ) +U IN ( 6 ) ) + 

# 24.*US35 + 28 . *U_IN (4) 

B = COF2 

COF2 = 9 . * ( -U IN ( 1 ) +U_IN ( 7 ) ) + 

# 6 . * ( -U IN ( 2 ) +U_IN ( 6 ) ) + 

# 3 . * ( -U_IN ( 3 ) +U_IN ( 5 ) ) 

q COF3 

COF3 = 5 . *US17 -3 . *US3 5 - 4.*U_IN(4) 

IF (ABS(COF3) .LT.1.0E-08) THEN 

PRTNT* ' ********* WARNING ********* 

PRINT*! 'CONSTANT VALUE EQUALS ZERO - CANNOT COMPUTE A 
PRINT* , ' FREQUENCY VALUE . ' 

RETURN 
END IF 

DID NOT DIVIDE BY 84 YET. DO SO FOR MAX PART BUT NOT PO. 

COMPUTE PO, VALUE OF F(P) WHERE A+B*P+C*P*P = 0 
COMPUTE FUNCTION AT PO; A+B*P0+C*P0*P0 = A-B**2/4C 
IF ( IOPT .EQ. 1) THEN 
PO=-0.5*COF2/COF3 

FMAX = (COFl + 0.5 * PO * COF2) /84. 


ELSE 
FMAX = 
END IF 
RETURN 
END 


(COFl + P0*( COF2 + PO*COF3 ) )/84 


ial = 7141, 
ia2 = 8121, 
ia3 = 4561, 


then 


icl = 54773, rml = 1 
ic2 = 28411, rm2 = 1 
ic3 = 51349) 


function ranl(idum) 
dimension r(97) 
parameter (ml = 259200, 

parameter (m2 = 134456, 

parameter (m3 = 243000, 

data iff /0/ 

if ( idum. It . 0 . or . if f • eq . 0) 
if f = 1 

ixl = mod (icl - idum, ml) 

ixl — mod(ial*ixl + icl, ml) 

ix2 = mod (ixl, m2) 

ixl = mod(ial*ixl + icl, ml) 

ix3 = mod (ixl, m3) 

do j = 1,97 

ixl = mod(ial*ixl + icl,ml) 

ix2 = mod(ia2*ix2 + ic2,m2) 

r ( j) = (float (ixl) + f loat (ix 2 ) *rm2) *rml 

end do 
idum = 1 
if 

= mod(ial*ixl + icl, ml) 


end 

ixl 


MAXIMUM' 


. 0/ml) 
. 0/m2) 
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ix2 = mod ( ia2*ix2 + ic2,m2) 
ix3 = mod(ia3*ix3 + ic3,m3) 
j = 1 + (97*ix3) /m3 
if ( j . gt . 97 . or . j . It . 1 ) pause 
rani = r(j) 

r ( j) = (float ( ixl) + f loat (ix2) *rm2) *rml 

return 

end 


32 



APPENDIX 2.E 


Simulation Test Results 


These are the results of simulation 3: 


MODEL DATA HAS 591 TIME POINTS 

MODEL DATA HAS 901 TIME POINTS 

WILL USE 591 TIME POINTS 

DT FOR MODEL DATA IS : 1.000000000 

DT FOR TRUTH DATA IS : 1.000000000 

OVERALL RMS MAGNITUDE ERROR = 34.42346191 

OVERALL RMS PHASE ERROR = 22.10563469 


The results of simulation 4 are as follows: 

MODEL DATA HAS 1601 TIME POINTS 

MODEL DATA HAS 1801 TIME POINTS 

WILL USE 1601 TIME POINTS 

DT FOR MODEL DATA IS : 1.000000000 

DT FOR TRUTH DATA IS : 1.000000000 

OVERALL RMS MAGNITUDE ERROR = 9.227395058 

OVERALL RMS PHASE ERROR = 6.498259544 


These are the results of simulation 5: 

MODEL DATA HAS 1601 TIME POINTS 

MODEL DATA HAS 1801 TIME POINTS 

WILL USE 1601 TIME POINTS 

DT FOR MODEL DATA IS : 1.000000000 

DT FOR TRUTH DATA IS : 1.000000000 

OVERALL RMS MAGNITUDE ERROR = 5.796232224 

OVERALL RMS PHASE ERROR = 4.697713375 


% 

DEGREES 


% 

DEGREES 


% 

DEGREES 



APPENDIX 2 . F 


Systematic^Tes^^ ^ No i se _f re e Test Cases, Station 2 

(2) NOISE - Noisy Test Cases, Station 2 

(3) LI BRAT I ON - Noise-free Tests at Station 1 

(4) LIB NOISE - Noisy Tests at Station 1 



The following 33 plots are the results of the 
NO NOISE. FOR program. (Refer to section I part C of the test 
plan.) These plots represent the amplitude, frequency, and phase 
errors for the eleven skiprope frequencies running from 0.0045 Hz 
to 0.0055 Hz vs. the penduluos and skiprope phases. 



Maximum Amplitude Error = 0.118% 
Frequency = 0.0045 Hz 




Error (per cent) xIO —1 
\O.Q O. 1 4 0.28 0.42 0.56 0.7 


Maximum Amplitude Error — 0.070% 
Frequency = 0.0049 Hz 



O'O 



Maximum Amplitude Error — 0.260% 
Frequency = 0.0051 Hz 




t) xIO 


Maximum Amplitude Error — 0.121% 
Frequency = 0.0052 Hz 



Maximum Amplitude Error = 0.266% 
Frequency = 0.0053 Hz 




cent 




Maximum Frequency Error — 0.434% 
Frequency = 0.0045 Hz 





Maximum Frequency Error — 0.590% 
Frequency = 0.0047 Hz 
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Maximum Frequency Error — 0.550% 
Frequency = 0.0053 Hz 




Maximum Frequency Error = 0.497% 
Frequency = 0.0055 Hz 







t) xIO -7 
9 1.2 1.5 


Maximum Phase Error — 0.147% 
Frequency = 0.0046 Hz 
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Maximum Phase Error = 0.127% 
Frequency = 0.0048 Hz 












